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ABSTRACT 

Recent observations support the hypothesis that a large fraction of "short-hard" gamma-ray bursts 
(SHBs) are associated with the inspiral and merger of compact binaries. Since gravitational-wave 
(GW) measurements of well-localized inspiraling binaries can measure absolute source distances with 
high accuracy, simultaneous observation of a binary's GWs and SHB would allow us to directly and 
independently determine both the binary's luminosity distance and its rcdshift. Such a "standard 
siren" (the GW analog of a standard candle) would provide an excellent probe of the relatively 
nearby (z < 0.3) universe's expansion, independent of the cosmological distance ladder, and thus 
complementing other standard candles. Previous work explored this idea using a simplified formalism 
to study measurement by advanced GW detector networks, incorporating a high signal-to- noise ratio 
limit to describe the probability distribution for measured parameters. In this paper we eliminate this 
simplification, constructing distributions with a Markov Chain Monte Carlo technique. We assume 
that each SHB observation gives both the source sky position and the time of coalescence, and we 
take both binary neutron stars and black hole-neutron star coalescences as plausible SHB progenitors. 
We examine how well parameters (particularly the luminosity distance) can be measured from GW 
observatations of these sources by a range of ground-based detector networks. We find that earlier 
estimates overstate how well distances can be measured, even at fairly large signal-to-noise ratio. 
The fundamental limitation to determining distance to these sources proves to be the gravitational 
waveform's degeneracy between luminosity distance and source inclination. Despite this, we find that 
excellent results can be achieved by measuring a large number of coalescing binaries, especially if 
the worldwide network consists of many widely separated detectors. Advanced GW detectors will be 
able to determine the absolute luminosity distance to an accuracy of 10-30% for NS-NS and NS-BH 
binaries out to 600 and 1400 Mpc, respectively. 

Subject headings: cosmology: distance scale — cosmology: theory — gamma rays: bursts — gravitational 
waves 



1. INTRODUCTION 
1.1. Overview 

There are presently two operational multikilometer in- 
terfcrometric gravitational-wave (GW) detectors: LIGO 4 
and Virgo 5 . They are sensitive to the GWs produced 
by the coalescence of two neutron stars to a distance of 
roughly 30 Mpc, and to the coalescence of a neutron star 
with a 10M Q black hole to roughly 60 Mpc. Over the 
next several years these detectors will undergo upgrades 
which are expected to extend their range by a factor 
~ 10. At these advanced sensitivity levels, most esti- 
mates suggest that detectors should measure at least a 
few, and possibl y a few dozen, binary co alescence events 
every year (e.g.. Kopparapu et al.ll2008l ). 

It has long been argued that neutron star-neutron star 
(NS-NS) and neutron star-black hole (NS-BH) merg- 
ers are likely to be accompanied by a gamma-ray burst 
(|Eichler et al.lfl989[ l. Recent evidence supports the hy- 
pothesis that many short-hard gamma-ra y bursts (SHBs) 
are indeed associated with such mergers jFox et al.ll2005t 
iNakar et al.ll2006l iBerger et al1l2007l . iPerlev et alJl2008ft . 
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This suggests the exciting possibility that it may be pos- 
sible to simultaneously measure a binary coalescence in 
gamma rays (and associated afterglow emission) and in 
GWs. The combined electromagnetic and gravitational 
view of these objects is likely to teach us substantially 
more than what we learn from either data channel alone. 
Because GWs track a system's global mass and energy 
dynamics, measuring GWs from a coalescing binary al- 
lows us to determine with exquisite accuracy "intrinsic" 
binary properties, such as the masses and spins of its 
members. As we describe in the following subsection, 
GWs can also determine a system's "extrinsic" prop- 
erties, such as location on the sky and distance to the 
source. In particular, the amplitude of a binary's GWs 
directly encodes its luminosity distance. Direct measure- 
ment of a coalescing binary could thus be used as a cos- 
mic distance measure: Binary inspiral would be a "stan- 
dard siren" (the GW equivalent of a standard candle, 
so-called due to the sound-like nature of GWs) whose 
calibr ation depends only on the validity of general rela- 
tivity (jDalal et al.ll2006l ). 

Unfortunately, GWs alone do not measure extrinsic 
parameters as accurately as the intrinsic ones. As we de- 
scribe in more detail in the following section, in general 
a GW observation of a binary measures a complicated 
combination of the distance to the binary, the binary's 
position on the sky, and the binary's orientation, with 
overall fractional accuracy ~ 1/signal-to- noise. As the 
distance is degenerate with the angular parameters, us- 
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ing GWs to measure absolute distance to a source re- 
quires a mechanism to break the degeneracy. Associat- 
ing the GW coalescence waves with a short-hard gamma- 
ray burst (SHB) is a near-perfect way to break this de- 
generacy. In this paper we explore the ability of near- 
future GW detector networks (such as LIGO and Virgo) 
to constrain binary parameters (and in particular, dis- 
tance to a binary), when used in conjunction with elec- 
tromagnetic observations of the same event (such as an 
associated SHB). We also examine how well these mea- 
surements can be improved if planned detectors in West- 
ern Australia (AIGO 6 ) and in Japan's Kamioka mine 
(LCGT 7 ) are operational. This paper substantially up- 
dates and improves upon earlier work palal et al.ll2006l . 
hereafter DHHJ06), utilizing a significantly more sophis- 
ticated and accurate parameter estimation technique. In 
the next section we review standard sirens, while in Sec. 
11.31 we briefly summarize DHHJ06. We follow this with 
a subsection describing the organization and background 
relevant for the rest of the paper. 

1.2. Standard sirens 

The GWs produced by the inspiral of two compact 
bodies directly encode the luminosity distance to the 
binary. It has long been recognized that GW in- 
spiral measurements could thus b e used as powerful 
tools for cosmology. ISchutzl (|1986l ) first demonstrated 
this by analyzing how binary coalescences a llow a di- 
rect m eas urement of the Hubble constant; iMarkovicI 
(|1993f) and iFinn fc Chernofil l)1993f) subsequently gener- 
alized this approach to include other cosmological pa- 
rameters. More recently, there has been much inter- 
est in the measurements enabled when the GWs from 
a merger are accompa nied by a counterpart in the elec- 
tromagnetic spectrum ([Bloom et afl feOOQ. Ph innevll2009l . 
iKulkarni fc Kasliwalll2009f) . In this paper we focus exclu- 
sively on GW observations of binaries that have an in- 
dependent sky position furnished by electromagnetic ob- 
servations (e.g., of gamma-rays associated with an SHB. 
or an accompanying optical afterglow). 

We begin by examining the gravitational waves from a 
binary inspiral as measured in a single GW detector. For 
simplicity, we only present here the lowest order contri- 
bution to the waves; in our subsequent calculations, our 
results are taken to higher order (see Sec. l2.1|) . The wave- 
form generated by a source at luminosity distance Dj,, 
corresponding to redshift z, is given by 
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Here <3>at is the orbital phase, / is the GW frequency, 
and M. = m\^ b rru/ b / {mi +m2) 1//5 is the binary's "chirp 
mass," which sets the rate at which 
changes. We use units with G = 1 = 



version factors are M, 
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Mq = GM Q /c 3 = 4.92 x 1CT 6 seconds. The angle i 
describes the inclination of the binary's orbital plane to 
our line-of-sight: cos i = L • n, where L is the unit vector 
normal to the binary's orbital plane, and n is the unit 
vector along the line-of-sight to the binary. The param- 
eters t c and <I> C are the time and orbital phase when the 
frequency / diverges in this model. (In reality, we expect 
finite size effects to substantially impact the waveform 
before this.) 

A given detector measures a linear combination of 
these polarizations: 

/Was = F+(9, 0, ^)h+ + F x (9, 0, iP)h x , (2) 

where 9 and </> describe the binary's position on the sky, 
and the "polarization angle" ip sets the inclination of 
the components of L orthogonal to n. The angles t and 
ip fully specify the orientation vector L. For a particu- 
lar detector geo metry, th e anten na functions F + and F x 
can be found in iThornel (|1987t) . (In Sec. H21 we give a 
general form for the gravitational waveform, without ap- 
pealing to a specific det ector, following the analysis of 
ICutler fc Flanagan! \l 994 hereafter abbreviated CF94.) 

Several features of Eqs. (TT|) and ^j) are worth com- 
menting upon. First, note that the phase depends on the 
redshifted chirp mass. Mea suring phase thus deter mines 
the redshifted chirp mass (jFinn fc ChernofflH99l . To 
understand why we measure redshifted chirp mass, note 
that M. controls how fast the frequency evolves: using 
Eq. we find / oc / n / 3 .M 5 / 3 . One can thus regard 
the chirp mass as entering the system's dynamics as a 
chirp time r c = GA4/c 3 . For a source at cosmological 
distance, this timescale is redshifted; the chirp mass we 
infer is likewise redshifted. Redshift and chirp mass are 
inextricably degenerate. Thi s remains true even when 
higher order effects (see, e.g.. iBlanchetj 120061 ) are taken 
into account: parameters describing a binary impact its 
evolutionary dynamics as timescales which suffer cosmo- 
logical redshift; the inferred values of those parameters 
are thus redshifted. GW observations on their own can- 
not directly determine a source's redshift. 

We next note that the amplitude depends on (l + z)M, 
the angles {9,<f>, and the luminosity distance D^. 

Measuring a GW amplitude thus measures some tan- 
gled combination of these parameters. By measuring 
the phase, we measure the redshifted chirp mass suffi- 
ciently well that (1 + z)A4 essentially decouples from 
the amplitude. More concretely, matched filtering the 
datastream with waveform templates should allow us to 
determine the phase with fractional accuracy ~ 
1/ [(signal-to- noise) x (number of measured cycles)]; (1 + 
z)M. should be measured with similar fractional accu- 
racy. Since NS-NS binaries will radiate roughly 10 4 cy- 
cles in band, and NS-BH binaries roughly 10 3 cycles, the 
accuracy with which phase and redshifted chirp mass can 
be determined should be exquisite. 

Although (1 + z)M. therefore decouples from the am- 
plitude, the distance, position, and orientation angles re- 
main highly coupled. If our goal is to determine the 
source distance, we must break the degeneracy that the 
amplitude's functional form sets among these parame- 
ters. For sources that we measure with ground-based 
GW detectors, one way to break these degeneracies is 
to measure the waves with mu l tiple detectors. Studies 
(|Svlvestrd l2TJ0l iCavalier et all [20061: iBlair et all 120081 ) 
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have shown that doing so allows us to determine the po- 
sition of a merging binary to within a few degrees, giving 
us some information about the source's distance and in- 
clination. 

Perhaps the best way to break these degeneracies 
would be to measure the event electromagnctically. An 
EM signature is almost certain to pin down the event's 
position far more accurately than GWs alone. The po- 
sition angles then effectively decouple, much as the red- 
shifted chirp mass decoupled. Using multiple detectors, 
we can then determine the source's orientation and its 
distance. This gives us a direct, calibration-free measure 
of the distance to a cosmic event. The EM signature 
may also provide us with the event's redshift, breaking 
the GW degeneracy of redshift and intrinsic parameters, 
and directly measuring both an event's distance and red- 
shift. In addition, if modeling or observations give us 
evidence for beaming of the SHB emission, that could 
strongly constrain the source inclination. 

1.3. This work and previous analysis 

Our goal is to assess how well we can determine the 
distance Dl to SHBs under the assumption that they 
are associated with inspiral GWs. We consider both 
NS-NS and NS-BH mergers as generators of SHBs, and 
consider several plausible advanced detector networks: 
the current LIGO /Virgo network, upgraded to advanced 
sensitivity; LIGO/ Virgo plus the proposed Australian 
AIGO; LIGO/Virgo plus the proposed Japanese LCGT; 
and LIGO/Virgo plus AIGO plus LCGT. 

The engine of our analysis is the construction of a 
probability function describing how measured source pa- 
rameters (focusing in particular on Dl) should be dis- 
tributed following GW measurement. Briefly, this for- 
malism builds a probability distribution for the parame- 
ters 9 of a source's waveform. (Components a of the vec- 
tor 9 are physical parameters such as a binary's masses, 
distance, sky position angles, etc.) Consider one detector 
which measures a datastream s(t), containing noise n(t) 
and a GW signal h(t,0), where 9 descr ibes the sou rce's 
"true" parameters. In the language of iFinnl (|1992t) . we 
assume "detection" has already occurred; our goal in this 
paper is to focus on the complementary problem of "mea- 
surement." Our aim is to assess how well we can deter- 
mine these par ameters from our data. 

As shown bv lFinnI (|1992f ) . given a model for our signal 
h(t, 9), and assuming that the noise statistics arc Gaus- 
sian, the probability that the parameters 9 describes the 
data s is 

p(9\ S )=po(9)cxp[-((h(9)-s)\(h(9)-s))/2} . (3) 
The inner product (a\b) describes the noise weighted 
cross-correlation of a(t) with b(t), and is defined precisely 
below. The distribution po(9) is a prior probability dis- 
tribution; it encapsulates what we know about our signal 
prior to measurement. We define 9 to be the parameters 
that maximize Eq. (|3j). 

DHHJ06 did a first pass on the analysis which we de- 
scribe here. They used an expansion in the variables 
(9 — 9) of the exponential to second order in the limit 
of large signal-to-noise ratio (SNR), which we will h ence- 
forth refer to as the "Gaussian" approximation (cf. IFinnl 
[1992): 

cxp[- (h(9) - s\h(9) - s) /2] ~ 



exp 




(4) 



where S8 a = 9 a - 8 a . In this limit, 9 = 9 (at least for 
uniform priors). The matrix 



(5) 



is the Fisher information matrix; 
covariance matrix. Diagonal entries E aa describe the 
variance of parameter 6 a ; off-diagonal entries describe 
correlations between different parameters. 

The Gaussian approximation to Eq. ([3]) is known to be 
accurate when the SNR is large; however, it is not clear 
what "large" really means. Given current binary coales- 
cence rate estimates, it is expected that most events will 
come from Dl ~ a few x 100 Mpc. In such cases, we can 
expect an advanced detector SNR ~ 10. It is unclear 
that this value is sufficiently high such that the "large 
SNR" approximation is appropriate. Indeed, some of the 
results in DHHJ06 seem somewhat anomalous. For ex- 
ample, it had been expected that an Australian detector 
would have a disproportionate impact on the network's 
ability to determine distance and inclination for many 
events, since a southern hemisphere detector's antenna 
functions would be substantially different than its north- 
ern cousins, offering excellent constraints on an event's 
position and inclination. In fact, DHHJ06 found that a 
southern detector's main impact was simply to add addi- 
tional SNR to a measurement — helpful, but not dispro- 
portionately so. We are concerned that this might be an 
artifact of the Gaussian approximation used in DHHJ06. 

In this analysis we do not use this Gaussian approx- 
imation, but instead use Markov Chain Monte Carlo 
(MCMC) techniques (in particular, the Metropolis- 
Hastings algorithm) to explore our parameter distribu- 
tions. The details of this techni que are summarized i n 
Sec. 3, and described in detail in lLewis fc Bridle! (|2002f K 
We find that MCMC techniques indeed indicate that the 
Gaussian approximation to Eq. ([3]) is failing in its esti- 
mate of a system's extrinsic parameters (though it ap- 
pears to do quite well for intrinisic parameters such as 
masses). In particular, we find that an Australian detec- 
tor would, in fact, have an outstanding impact on our 
ability to use SHB GWs as standard sirens. In addition, 
DHHJ06 restrict themselves to NS-NS binaries, observed 
with a four detector network. In what follows we relax 
both of these assumptions. 

1.4. Organization of this paper 

We begin in Sec. [5] by summarizing how GWs encode 
the distance to a coalescing binary. We first describe 
the post-Newtonian (PN) gravitational waveform we use 
in Sec. 12.11 and then describe how that wave interacts 
with and can be measured by a network of detectors in 
Sec. 12.21 Our discussion of the network-wave interaction 
is heavily based on the notation and fo rmalism used in 
Sec. 4 of CF94, as well as the analysis of lAnderson et ahl 
(|2001h . Section 12.21 is sufficiently dense that we sum- 
marize its major points in Sec. 12.31 before concluding, in 
Sec. 12.41 with a description of the various GW detectors 
which we include in our analysis. 
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We outline our parameter estimation formalism in Sec. 
3. In Scc. l3.ll we describe in more detail how to construct 
the probability distributions describing parameter mea- 
surement mentioned above (Sec. II. 3[) . We then give, in 
Sec. 13.21 a brief description of our selection procedure 
based on SNR detection thresholds. This procedure sets 
physically motivated priors for some of our parameters. 
The Markov Chain Monte Carlo technique we use to ex- 
plore this function is described in Sec. 13.31 How one ap- 
propriately averages this distribution to give "noise aver- 
aged" results (and to compare with previous literature) 
is discussed in Sec. 13.41 

In Sec. H] we discuss the validation of our code. We be- 
gin by attempting to reproduce some of the key results 
on distance measurement presented in CF94. Because of 
the rather different techniques used by Cutler & Flana- 
gan, we do not expect exact agreement. It is reassuring 
to find, nonetheless, that we can reconstruct with very 
good accuracy all of the major features of their analysis. 
We then examine how these results change as we vary 
the amplitude (moving a fiducial test binary to smaller 
and larger distances), as we vary the number of detectors 
in our network, and as we vary the source's inclination. 

Our main results are given in Sec. [5] We consider sev- 
eral different plausible detector networks and examine 
measurement errors for two "fiducial" binary systems, 
comprising either two neutron stars (NS-NS) with phys- 
ical masses (i.e., rest frame masses, not including red- 
shift effects) of mi = mi = 1.4 Mq, or a neutron star 
and black hole (NS-BH) system with physical masses 
mi = 1.4 Mq and 777,2 = 10 M@. Assuming a constant 
comoving cosmological density, we distribute potential 
GW-SHB events on the sky, and select from this distribu- 
tion using a detection threshold criteria set for the entire 
GW detector network. We summarize some implications 
of our results in Sec. [H A more in-depth discussion of 
these implications, particularly with regard to what they 
imply for the ability of standard sirens to measure the ex- 
pansion of the universe, will be presented in a companion 
paper. 

Throughout this paper we use units with G = 1, c = 1. 
Because we often refer to "redshifted" masses, we define 
the shorthand m z = (1 + z)m for any mass parameter 
777. 

2. MEASURING GRAVITATIONAL WAVES FROM 
INSPIRALING BINARIES 

In this section we review the GW description used in 
our analysis, the formalism describing how these waves 
interact with a network of detectors, and the properties 
of the detectors. 

2.1. GWs from inspiraling binaries 

The inspiral and merger of a compact binary's mem- 
bers can be divided into three consecutive phases. The 
first and longest is a gradual adiabatic inspiral, when 
the members slowly spiral towards one another driven 
by the loss of orbital energy and angular momentum due 
to radiative backreaction. Post-Newtonian (PN) tech- 
niques (an expansion in the gravitational potential M/r, 
or equivalently for bound systems, the orbital speed v 2 ) 
allow a binary's evolution, and its emi tted GWs , to be 
modeled analytically to high order; see iBlanchetl (|2006f ) 
for a review. When the bodies come close together, the 



PN expansion is no longer valid, and direct numerical 
calculation is required. Recent breakthroughs in numeri- 
cal relativity are making it possible to fully model the 
strong-fi eld, dynamical merger of the two bodie s into 
one; see IPretoriusI (12001 . IShibata fc Ur^ol (l200l . and 
lEtienne et al.l (|2008f ) for discussion. If the end state is a 
single black hole, the final waves from the system should 
be described by a ringdown as the black hole settles down 
to the Kerr spacetime solution. Note that much success 
has been achieved in analytically modeling the entire coa- 
lescence process usi ng the "effective one body" tec hnique. 
Initially proposed in lBuonanno fc Damouri (|1999l ). recent 
results show good agreement over the entire coalescence 
wit h the most accurately a vailable numerical simulations 
(see IBuonanno et al.ll2009T) . 

In this work we are concerned only with the inspiral 
phase, and will accordingly use the PN waveform to de- 
scribe our waves. In particular, we use the so-called "re- 
stricted" PN waveform. Following CF94, the waveform 
of a binary inspiral may be written schematically 




(0) 



Here x indicates the PN order [h x is computed to 0{v 2x ) 
in orbital speed], m denotes harmonic order (e.g., 777 = 2 

is quadrupole), and $ or h{t) = j £l(t')dt' is orbital phase 
[with f2(<) the orbital angular frequency]. The "re- 
stricted" waveform neglects all PN amplitude terms be- 
yond the leading one, and considers only the dominant 
77i = 2 contribution to the phase. The phase is itself 
computed to high PN order. As argued in CF94, these 
restrictions are useful because accumulated phase infor- 
mation has a greater impact on measurement accuracies 
than information in the amplitude. 

Let the unit vector n point to a binary on the sky 
(so that the waves propagate to us along — n), and let 

the unit vector L denote the normal along the binary's 
orbital angular momentum. The waveform is fully de- 
scribed by the two polarizations: 

h+(t) = ^ {irM z f(t)} 2/3 [1 + (L • n) 2 ] cos[<D(t)] , 



hy.it) = 



D L 

D L 
4M 



{irM z f(t)f 3 _4+(n,L)cos[$(i)] ; (7) 
{TrM z f(t)}V 3 (L-h)sm{<f>(t)} , 



[7rM z f(t)f 3 _A x (n,L)sin[$(i)] . (8) 

J->L 

Equations ([7]) and ([5]) are nearly identical to those given 
in Eq. {T]). In particular, M. z is the binary's redshifted 
chirp mass, Dl is the luminosity distance to the binary, 
and we have explicitly defined the inclination angle cos 1 
in terms of the vectors n and L. For later convenience, 
we have moved all dependence on sky position and ori- 
entation into the functions A+. x - In Sec. 12.21 we discuss 
how these polarizations enter the GW field more gener- 
ally, and how they interact with our detectors. 

A major difference in these forms of h + and h x , as 
compared to Eq. (fT]), is that we now compute the GW 
phase to higher order: 
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using the 2nd-post-N ewtonian (2PN) frequency sweep 
(jBlanchet et alJll995[ ) 

dt 5 z J 
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Higher order results for d f / dt are now known 
(jBlanchet et alJ l2002al lE l2004h . but 2PN order will be 
adequate for our purposes. Since distance measurements 
depend on accurate amplitude determination, we do not 
need a highly refined model of the wave's phase. The rate 
of sweep is dominantly determined by the chirp mass, 
but there is an important correction due to 77 = [i/M = 
miiri2/(mi + 7712) 2 , the reduced mass ratio. Note that 77 
is not redshiftcd; both /j, and M (the reduced mass and 
total mass, respectively) acquire (1 + z) corrections, so 
their ratio is the same at all z. Accurate measurement of 
the frequency sweep can thus determine both A4 Z and 77 
(or M z and /i z ), potentially determining the redshiftcd 
masses of each member of the binary. 

We will often find it useful to work in the frequency 
domain, using the Fourier transform h(f) rather than 
hit): 



Hf) 



,2mft 



h(t) dt 



(11) 



An approximate analytical result for h(f) can be 
found using the sta tionary phase approximation 
(|Finn fc Chernofj |1993[ ). which describes the Fourier 
transform when / changes slowly: 
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"Slowly" means that / does not change very much over 
a single wave period 1//, so that (df /dt)/ f <C /. The 
validity of this approximation for the waveforms we con- 
sider, at least until the last moments before merger, ha s 
been demonstrated in previous work (|Droz "eFalJ[l999h . 
The phase function *(/) in Eqs. (jT2j) and ([13]) is given 
by 
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As in Eq. (JTJ) , t c is called the "time of coalescence" and 
defines the time at which / diverges within the PN frame- 
work; <E> C is similarly the "phase at coalescence." Follow- 
ing the method employed in earlier works, we assume an 
abrupt (and unphysical) transition between the inspiral 
and merger phases at the so-called innermost stable cir- 
cular orbit (ISCO), /isco = (GV&kMz)- 1 . For NS-NS, 



/isco occurs at high frequencies where ground-based de- 
tectors have poor sensitivity; as such, we are confident 
that this abrupt transition has little impact on our re- 
sults. For NS-BH, /isco is likely to be in a band with 
good sensitivity. Better modeling of this transition could 
be important in this case. 

In this analysis, we neglect effects which depend on 
spin. In general relativity, spin drives precession effects 
which can "color" the waveform in important ways, and 
which can in principle have important observational ef- 
fects (see, e.g.. lLang fc H ughes 2006, van der Sluvs et aTl 
l2008h . These effects are important when the dimension- 
less Kerr spin parameter is fairly large. Neutron stars are 
unlikely to spin fast enough for their angular momentum 
to drive interesting precessions during the time that they 
are in the band of GW detectors. To show this, write the 
moment of inertia of a neutron star as 



NS-n-NS > 



(15) 



where Mns and i?NS are the star's mass and radius, and 
the parameter k describes the extent to which its mass 
is centrally condensed (compared to a uniform sphere). 
Detailed calculation s with differ e nt equ ations of state in- 
dicate k ~ 0.7-1 [cf. lCook et alJ (j!994f ): this range comes 
the slowly rotating configurations in their Tables 12, 15, 
18, and 21]. For a neutron star whose spin period is Pns, 
one finds that the Kerr parameter is given by 



2ncL 



NS 



GM* s Pss 



:0.06k 



-Rns 
12 km 



1.4 M, 



xs 



10 msec 
Pt 



xs 



(16) 



As long as the neutron star spin period is longer than 
~ 10 msec, <2ns is small enough that spin effects can 
probably be safely neglected in our analysis. Strictly 
speaking, we should include spin in our models of BH- 
NS binaries; we leave th i s to a later analysis. We note 
that Ivan der Sluvs et"ail (|2008f ) included black hole spin 
effects in an analysis which did not assume known source 
position. They found that spin-induced modulations 
greatly improved the ability of GW detectors to local- 
ize a source, even when measured by only one detector. 
This suggests that, if position is known, spin modula- 
tions would likewise greatly improve our ability to mea- 
sure source inclination and distance. 

Our GWs depend on nine parameters: two mass pa- 
rameters M. z and fx z , two sky position angles (which 
set n), two orientation angles (which set L), time at 
coalescence t c , phase at coalescence $ c , and luminos- 
ity distance to the source D^. When sky position is 
known (e.g., from an SHB observation, as assumed in 
this analysis), the parameter set is reduced to seven: 
{M z , \l z , D L , t c , cos 1, tp, $ c }. 

2.2. Measurement of GWs by a detector network 

We now examine how the waves described in Sec. I2.ll 
interact with a network of detectors. We begin by intro- 
ducing some conventions to describe our binary and our 
detectors. We do this geometrically, using vectors which 
describe our source's position and orientation, and our 
detectors' positions and orientations. 
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As described in the previous section, a source's sky 
position is given by a unit vector n (which points from the 
center of the Earth to the binary), and its orientation is 
given by a unit vector L (which points along the binary's 
orbital angular momentum). We construct from these a 
pair of axes which describe the binary's orbital plane: 



n x L 
In x LI 



Y 



n x X 
In x X| 



(17) 



These axes in turn are used to define the polarization 
basis tensors 



e+=X(g)X-Y(g)Y 
e x = ±(g> Y + Y(X>X 



(18) 
(19) 

From these, the transverse-traceless metric perturbation 
which describes the GWs emitted by our source can be 
written 



h + e 



(20) 



We next characterize the GW detectors. Each detector 
is an L-shaped interferometer whose arms define two- 
thirds of an orthonormal triple. Denote by x a and y a 
the unit vectors along the arms of the a-th detector in 
our network; we call these the x- and y-arms. (The vec- 
tor z a = x a x y a points radially from the center of the 
Earth to the detector's vertex.) These vectors define the 
response tensor for detector a: 

Dj? = \ [(x a r(x Q r - (y a ) W] . (21) 
The response of detector a to a GW is given by 



K = D i Jh ij 



-27ri(n-x a )/ 



(F a .+h + + F a , x h x ) 



(22) 



where x a is the position of the ath detector and the factor 
(n ■ x a ) measures the time of flight between detector a 
and the coordinate origin. The second form of Eq. (|22|) 
shows how the antenna response functions introduced in 
Eq. ^ are built from the wave tensor and the response 
tensor. 

Thus far this discussion has been completely frame- 
independent, in that we have defined all of our vectors 
and tensors without reference to a particular coordinate 
system. To facilitate an analysis of how inspiral GWs 
interact with a network of detectors, we now introduce a 
co ordinate system f or our detectors following the analysis 
of I Anderson et all (120011) [who in turn use the WGS-84 
Earth model ( Althouse et al.|[l998f )]. In this system, the 
Earth is taken to be an oblate ellipsoid with semi-major 
axis a = 6.378137 x 10 6 meters, and semi-mmor axis 
b = 6.356752314 x 10 6 meters. Our coordinates are fixed 
relative to the center of the Earth. The x-axis (which 
points along i) pierces the Earth at latitude 0° North, 
longitude 0° East (normal to the equator at the prime 
meridian); the y-axis (along j) pierces the Earth at 0° 
North, 90° East (normal to the equator in the Indian 
ocean somewhat west of Indonesia) ; and the z-axis (along 
k) pierces the Earth at 90° North (the North geographic 
pole). 

A GW source at (9, <p) on the celestial sphere has sky 
position vector n: 



The polarization angle, ip, is the angle (measured clock- 
wise about n) from the orbit's line of nodes to the 
source's X-axis. In terms of these angles, the vectors 
X and Y are given by ((Anderson et al.ll2001h 

X = (sin <fi cos ip — sin ip cos </> cos 9)i 

— (cos <p cos ip + sin ip sin <p cos 9)j + sin ip sin 0k , 

(24) 

Y = (— sin cp sin ip — cos ip cos (p cos 9)i 

+(cos (p sin ip — cos ip sin <p cos 9)j + cos ip sin 9k . 

(25) 

Note that the angle (p is related to the right ascension, a, 
by a = </> + GMST (where GMST is the Greenwich mean 
sidereal time at which the signal arrives), and 9 is related 
to th e declination, <5, by 5 — n/2 — 9 (cf. I Anderson et alJ 
120011 Appendix B). Combining Eqs. ((24j) and ([25]) with 
Eqs. (fT^)) - (|2"0|) allows us to write hij for a source in co- 
ordinates adapted to this problem. 

We now similarly describe our detectors in terms of 
convenient coordinates. Detector a is at East longitude 
A a and North latitude (p a (not to be confused with sky 
position angle <p). The unit vectors pointing East, North, 
and Up for this detector are 

e^ = -sinA a i + cosA a j , (26) 

e a = — sin ip a cos A a i - sin ip a sin A a j + cos (p a k , (27) 

= cos tp a cos A a i + cos ip a sin A a j — cos </? a k . (28) 

The x-arm of detector a is oriented at angle T a North 
of East, while its y-arm is at angle T a + 7r/2. Thanks to 
the Earth's oblateness, the x- and y-arms are tilted at 
angles u)%' v to the vertical. The unit vectors x a , y a can 
thus be written 

x a =cosw a cosT a e a +cosw a sinT a e a +smu a e , 

(29) 

y a = - cos uj v a sin T a ef + cos uj y a cos T a e^ + sin u^e u . 

(30) 

Combining Eqs. (J29J) and (J30J) with Eq. (J2TJ) allows us to 
write the response tensor for each detector in our net- 
work. 

2.3. Summary of the preceding section 

Section 12.21 is sufficiently dense that a brief summary 
may help clarify its key features, particularly with re- 
spect to the quantities we hope to measure. From Eq. 
(|22[) , we find that each detector in our network measures 
a particular weighted sum of the two GW polarizations 
h + and h x . We can rewrite the waveform detector a 
measures as 

h a = {irM z f(t)f 3 cos [$(i) + %} , (31) 

where we have introduced detector a's "polarization am- 
plitude" 



Ap= {F^+A+f + (F a , x A x ) 2 , 
and its "polarization phase" 

F a ,xA x 



n = sin 9 cos (pi + sin 9 sin <pj + cos 6>k 



(23) 



tan $ T 



Fa, + A+ 



(32) 



(33) 
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I Cut led (|1998f ) introduced these quantities to describe the 
response of the space-based LISA detector to a source's 
waves. As previously discussed, the intrinsic GW phase, 
is a strong function of the redshifted chirp mass, 
M. Zl the reduced mass, fj, z , the time of coalescence, t c , 
and the phase at coalescence, $ c . Measuring the phase 
determines these four quantities, typically with very good 
accuracy. 

Consider for a moment measurements by a single de- 
tector. The polarization amplitude and phase depend 
on the binary's sky position, (0, cp) or n, and orienta- 
tion, [ip, l) or L. [They also depend on detector position, 
(\ a ,(p a ), orientation, T a , and tilt, (w£,cj^). Since these 
angles are known and fixed, we ignore them in this dis- 
cussion.] If the angles (6,<p,ip,b) are not known, then a 
single detector cannot separate them, nor can it separate 
out the distance, Dl- 

Multiple detectors can, at least in principle, separately 
determine these parameters. Each detector measures its 
own amplitude and polarization phase, and by combin- 
ing their outputs, we can fit to the four unknown an- 
gles and the distance, Dl- Various works have analyzed 
how well this can be done assuming th at the position 
and orientation are completely unknow n (jSvlvestrdl2004l : 
ICavalier et alll2006t UTair et al.ll2008f) . Van der Sluys et 
al. (2008) performed such an analysis for measurements 
of NS-BH binaries, including the effect of orbital preces- 
sion induced by the black hole. These precessions effec- 
tively make the angles u and tp time dependent, breaking 
the degeneracy among these angles and Di in a com- 
pletely different way. 

In what follows, we assume that an electromagnetic 
identification pins down the angles (0,4>), so that they 
do not need to be determined from the GW data. We 
then face the substantially less challenging problem of 
determining ip, l, and Dl from the various waveforms, 
h ai that we measure at the detectors. We will also exam- 
ine the impact of a constraint on the inclination, l. Long 
bursts are believed to be strongly collimated, emitting 
into jets with opening angles of just a few degrees. Less 
is known about the collimation of SHBs, but it is plausi- 
ble that their emission may be primarily along a preferred 
axis (presumably the progenitor binary's orbital angular 
momentum axis). We are particularly interested in as- 
sessing the extent to which constraining t to be within 
some Gaussian distribution with opening angle 8l ~ 25° 
of the "face-on" configuration (|n • L| = 1) improves our 
ability to determine the source's distance. 

2.4. GW detectors used in our analysis 

Here we briefly summarize the properties of the GW 
detectors that we consider: 

LIGO: The Laser Interferometer Gravitational-wave Ob- 
servatory consists of two 4 kilometer interferometers 
located in Hanford, Washington (US) and Livingston, 
Louisiana (US). These instruments have achieved their 
"initial" sensitivity goals; an upgrade to "advanced" 
sensitivity is expected to be completed around 2015. 
We show the projected advanced sensitivity in Fig. [TJ 
Broadly speaking, the noise spectrum is dominated by 
seismic vibrations due to ground and human activity at 
/ < 10 Hz, thermal vibrations in the band 10 Hz < / < 
50 Hz, and photon shot noise for / > 50 Hz. At present, 




1 10 100 1000 10000 

frequency [Hz] 

Fig. 1. — Projected noise curve for Advanced LIGO. We assume 
a lower frequency cut-off of 10 Hz. 

the Hanford site also houses a 2 kilometer interferome- 
ter; this will be extended to a full 4 kilometer baseline 
in the advanced era. As such, the LIGO network will 
consist of two instruments at Hanford and one in Liv- 
ingston. Because the two Hanford instruments arc likely 
to have somewhat correlated noise (particularly at low 
frequencies), we assume only one instrument in Hanford 
for most of our analysis. 

Virgo: The Virgo detector near Pisa, Italy has slightly 
shorter arms than LIGO (3 kilometers), but should 
achieve similar advanced detector sensitivity, on roughly 
the same timescale as the LIGO detectors. For simplic- 
ity, we will assume that Virgo's sensitivity is the same as 
LIGO's in our analysis. 

Our baseline detector network consists of the LIGO 
Hanford and Livingston sites, and Virgo; these are in- 
struments which are running today, and will be upgraded 
over the next decade. We also examine the impact of 
adding two proposed interferometers to this network: 
AIGO: The Australian International Gravitational Ob- 
servatory is a proposed multikilometer interferometer 
that would be located in Gingin, Western Australia. 
AIGO's proposed site in Western Australia is particu- 
larly favourable due to low seismic and human activity. 
LCGT: The Large-scale Cryogenic Gravitational-wave 
Telescope is a proposed multikilometer interferometer 
that would be located in the Kamioka observatory, 1 kilo- 
meter underground. This location takes advantage of the 
fact that local ground motions tend to decay rapidly as 
we move away from the Earth's surface. They also plan 
to use cryogenic cooling to reduce thermal noise. 

As with Virgo, we will take the sensitivity of AIGO 
and LCGT to be the same as LIGO for our analysis. 
Table [1] gives the location and orientation of these detec- 
tors, needed to compute the response function for each 
member of our network. 

3. ESTIMATION OF BINARY PARAMETERS 
3.1. Overview of formalism 

We now give a brief summary of parameter estima- 
tion. In particular, we lay the foundations for estimat- 
ing measurement accuracies given a GW de tector ' s nois e 
spectrum. Further details are discussed in iFinnl (|1992f) . 
iKrolak et all (fl993l ). and CF94. 

The datastream of detector a, s a (t), has two major 
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TABLE 1 

GW DETECTORS WE CONSIDER (POSITIONS AND ORIENTATIONS) 



Detector 


East Long. A 


North Lat. ip 


Orientation T 


X 


arm tilt uj x 


j/-arm tilt u) v 


LIGO-Han 


-119.4° 


46.5° 


126° 


(- 


-6.20 X 10~ 4 )° 


(1.25 X 10- 5 )° 


LIGO-Liv 


-90.8° 


30.6° 


198° 


(- 


-3.12 X 10" 4 )° 


(-6.11 X 10" 4 )° 


Virgo 


10.5° 


43.6° 


70° 


0.0° 


0.0° 


AIGO 


115.7° 


-31.4° 


0° 




0.0° 


0.0° 


LCGT 


137.3° 


36.4° 


25° 




0.0° 


0.0° 



contributions: The true GW signal h a (t; 9) (constructed 
by contracting a source's GW tensor hij with detector 
a's response tensor D l J; cf. Sec. 12. 2[) . and a specific real- 
ization of detector noise n a (t), 

s a (t) = h a (t;0)+n a (t) . (34) 

The incident gravitational wave strain is assumed to de- 
pend on (unknown) true parameters 9. As in Sec. 11.31 
9 is a vector whose components are binary parameters. 
Below we use a vector s whose components s a are the 
datastreams of each detector (and likewise h and n are 
vectors whose components are the GW and noise content 
of each detector, respectively). 

We assume the noise to be stationary, zero mean, and 
Gaussian. This lets us categorize it using the spectral 
density, as follows. First, define the noise correlation 
matrix: 

= (n a (t + T~)n b (t)) - (n a (t + r)) (n b (t)} 

= (n a {t + T)n b (t)) , (35) 

where the angle brackets are ensemble averages over noise 
realizations, and the zero mean assumption gives us the 
simplified form on the second line. For a = b, this is 
the auto-correlation of detector a's noise; otherwise, it 
describes the correlation between detectors a and b. The 
(one-sided) power spectral density matrix is the Fourier 
transform of this: 



/oo 
dre 2 ^ fT C n (T) ab 
-OO 



(36) 



Note that this is defined for / > only. For a = b, this 
is the spectral density of noise power in detector a; for 
a 7^ b, it again describes correlations between detectors. 
From these definitions, one can show that 



(n a (f)Mf 



= \s( f 



f')S n {f)ab- 



(37) 



For Gaussian noise, this statistic completely character- 
izes our detector noise. No real detector is of course 
completely Gaussian, but by using multiple, widely- 
separated detectors non-Gaussian events can be rejected. 
In our analysis, we assume the detectors' noises are un- 
correlated such that Eq. (|37|) becomes 



(38) 



(M/) MfT) = \s ab s(f - f)s n (f) a - 



Finally, for simplicity we assume that S n (f) a has the uni- 
versal shape S n (f) projected for advanced LIGO. This 
spectrum is shown in Fig. [T] 

The central quantity of interest in parameter estima- 
tion is the posterior probability distribution function 
(PDF) for 9 given detector output s, which is defined 
as 

P (9\s) =N P {0 \9)£ TO t(s\9). (39) 



(Note that we assume a GW has been detected.) In 
Eq. ([39)) . TV denotes a normalization constant, p i -°' , (9) 
is the PDF that represents the prior probability that 
a measured GW is described by the paramete rs 9, and 
£tot (s I 9) is the total likelihood function (e.g.. lMacKavl 
120031 ). The likelihood function measures the relative con- 
ditional probability of observing a particular set of data 
s given that a measured signal h depending on some un- 
known set of parameters, 9 and noise, n. Because we 
assume that the noise is independent and uncorrelated 
at each detector site, we may take the total likelihood 
function to be the product of the individual likelihood at 
each detector: 



£tot(s|0) = U a C a (s a \9) 



(40) 



wh ere C a , the likelihood computed at detector a, is given 
by (|Finnlll992D 



C a (S I 9) = 

The inner product (. . 
is defined as 



(g\h) 



df 



,-(h a (e)-8 a \h a (e)- Sa )/2 _ / 41 x 
I . . .) on the vector space of signals 

~g*(f)kf) + ~g(f)h*(f) (A9] 

Sn(f) ■ { ' 

This definition means that the probability of the random 
noise n(t) assuming a particular realization no(i) is 

p(n = n ) cx e -(«ol"o)/2^ (43) 

For clarity, we distinguish between various definitions 
of the signal-to-noise ratio (SNR). The true SNR at de- 
tector a, associated with a given instance of noise for a 
measurement at a particular detector, is defined as (see 
CF94) 

S\ (h a \s a ) 

. N J a, true V( h " I h a) 

This is a random variable with Gaussian PDF of unit 
variance. For an ensemble of realizations of the detector 
noise n a , the so-called average SNR at detector a is given 

by 

S\ = (hg\K) 

N 1 rms (h a \n a ) 



(44) 



(fc«|fca) 1/2 . (45) 



Consequently, we can define the combined true and av- 
erage SNRs of a coherent network of detectors: 



and 



(46) 



(47) 
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Estimating the parameter set © often uses the "max- 
imum likelihood" method followi ng either a Bayesian 
(lLoredol U989L iFirml [liM C F94. iPoisson fe Willi Tl995h 
or frequcntist point of view (|Kr61ak et al.lll993l CF94) . 
We do not attempt to review these differing philosophies, 
and instead refer to Appendix A2 of CF94 for detailed 
discussion. It is worth noting that, in the GW literature, 
the "maximum likelihood" or "maximum a posterior" are 
often interchangeably referred to as the "best-fit" param- 
eters, which can be confusing. The maximum a poste- 
rior is the parameter set ©map which maximizes the full 
posterior probability, Eq. (|39p : likewise, the maximum 
likelihood is the parameter set ©ml which maximizes the 
likelihood function, Eq. ((40]) . 

Following the approach advocated by CF94, we intro- 
duce the Bayes estimator #bayes( s )i defined as 

£bayes(s) = J e l p(6\s)de. (48) 

The integral is performed over the whole parameter set 
©; dO = dd 1 dQ 2 . . . d6 n . Similarly, we define the rms 
measurement errors ^bayes 

^AYES = / P'-flWs) (^-^BAYEsMe I (49) 

The meaning of the Bayes estimator can be understood 
as follows. Consider a single detector which records an 
arbitrarily large ensemble of GW signals. Within this en- 
semble, there will be a subset of signals with the same de- 
tector output s(t). Although these outputs are identical, 
each actually corresponds to GW signals with different 
true parameters ©; it just so happens that they also have 
different noise realizations n(t) that conspire so that the 
sum s(t) is the same. In this case, #bayes( s ) re P resen ts 
the expectation of l averaged over the ensemble of GW 
signals. The principle disadvantage of the Bayes estima- 
tor is computation cost, owing to the multi-dimensional 
integrals in Eqs. (|4"8"]l and (|49|) . In Sec. 13.31 we describe 
MCMC methods which efficiently approximate these in- 
tegrals. 

For large SNR it can be shown that the estimators 
©ML; ©MAP: an d ©bayes agree with one another (CF94). 
Furthermore, in this limit (and as mentioned in Scc. ll.3[) 
Eq. (|39p assumes a Gaussian form with respect to param- 
eter errors, allowing us to estimate parameter variances 
as the inverse of a Fisher matrix. However, as illustrated 
in Sec. IVD of CF94, effects which scale nonlinearly with 
1/SNR and prior information [represented by the distri- 
bution p(°>(6)] contribute significantly at low SNR. The 
Gaussian approximation to Eq. (|39p then tends to un- 
derestimate measurement errors by missing tails or mul- 
timodal structure in posterior distributions. 

We emphasize that in this analysis we do not con- 
sider systematic errors that occur both due to limita- 
tions in our source model and due to gravitational lcns- 
ing effects. A framework for analyzing systematic er- 
ror s in GW measurem e nts h as recently been presented 
by iCutler fc Vallisneril (|2007| ). An important follow-on 
to this work will be to estimate systematic effects and 
determine whether they significantly change our conclu- 
sions. 



3.2. Binary Selection and Priors 

We now detail our method for generating a sample of 
detectable GW-SHB events. The selection is central to 
the use of GW-SHBs as standard sirens. In addition, 
it sets physically motivated priors for the parameters in 
our Bayesian framework. We employ a simple selection 
procedure based on detection thresholds. An implicit 
assumption of our method, detailed below, is that every 
SHB event is a GW candidate potentially detectable by 
a network of advanced GW interferometers; 

We assum e a constant comoving density (|Peeblesj|1993l 
lHog3[l999l ) of GW-SHB events, in a ACDM Universe 
with H n = 70.5 km/sec/Mp c, J7a = 0.726, and fl m = 
0.2732 (jKomatsu et al.|[2009h ; our results are insensitive 
to the precise details of the cosmology. We Monte-Carlo 
distribute 10 6 binaries uniformly in volume, with ran- 
dom sky positions and orientations, to redshift z = 1 
(i.e., to roughly 6.6 Gpc). We compute the average SNR, 
Eq. (|4~5|) . for each binary at each detector, and use Eq. 
(|47|) to compute the average total SNR for each network 
we consider. Since we assume prior knowledge of the 
merger time (by assumption that the inspiral is corre- 
lated with a SHB), we set a threshold SNR for the to- 
tal detector network, SNR to tai = 7.5 (see discussion in 
DHHJ06). This is somewhat reduced from the thresh- 
old we would set in the absence of a counterpart, since 
prior knowledge of merger time and source position re- 
duces t he number of search temp l ates we need by a factor 
~ 10 5 (jKochanek fc Piradll993L IOwenlll996f ). By using 
the average SNR to set our threshold, we introduce a 
slight error into our analysis, since the true SNR will dif- 
fer from the average. Some events which we identify as 
above threshold could be moved below threshold due to 
a measurement's particular noise realization. However, 
some sub-threshold events will likewise be moved above 
threshold due to the instance of noise. The net effect is 
not expected to be significant, and will be evaluated in 
future work. 

Our choice of threshold selects detectable GW-SHB 
events for each detector network. To avoid confusion, 
we use the term "total detected binaries" to mean bina- 
ries which are detected by a network consisting of all five 
detectors— both LIGO sites, Virgo, AIGO, and LCGT. 
As one might imagine, this five detector network detects 
more binaries than one with four detectors, which in 
turn detects more than one of three detectors. Both the 
Southern Hemisphere detector AIGO and the Japanese 
detector LCGT substantially improve on the number of 
binaries detected, as compared to the two LIGO detec- 
tors plus Virgo. Assuming isotropy in SHB orientation 
(i.e., that all binary orientations are equally likely given 
an SHB), we find that a LIGO-LIGO- Virgo network de- 
tects 50% of the total detected binaries; LIGO-LIGO- 
Virgo-AIGO detects 74% of the total; and LIGO-LIGO- 
Virgo-LCGT detects 72% of the total. Figure [2] shows 
the location of detected binaries on the sky for LIGO 
and Virgo (upper left panel); LIGO, Virgo, and AIGO 
(upper right); LIGO, Virgo, and LCGT (lower left); and 
all detectors (lower right). Notice that networks which 
include LCGT tend to have rather uniform sky cover- 
age, and that with AIGO, the quadrants with cos6> > 0, 
4> > it and with cos 9 < 0, <p < tt are covered particularly 
well. 
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LIGO+Virgo 



LIGO+Virgo+AIGO 




Fig. 2. — Detected NS-NS binaries for our various detector net- 
works, as a function of sky position (cos 6,0). The lower right 
panel shows all of the binaries detected by a five-detector network 
consisting of LIGO Hanford, LIGO Livingston, Virgo, AIGO, and 
LCGT; upper left is LIGO plus Virgo, with 50% of the five-detector 
events detected; upper right is LIGO, Virgo, and AIGO, with 74% 
of events detected; and lower left is LIGO, Virgo, and LCGT, with 
72% of events detected. Notice that detected binaries are more uni- 
formly distributed on the sky in networks that include LCGT, and 
that AIGO markedly improves the coverage in two of the sky's 
quadrants. Our coordinate is related to right ascension a by 
= o-GMST, where GMST is Greenwich Mean Sidereal Time; 
is related to declination (5by0 = 7r/2 — S. 

Our selection method implicitly sets a prior distribu- 
tion on our parameters. For example, the thresholding 
procedure results in a significant bias in detected events 
toward face-on binaries, with L-n — > ±1. Figure [3] shows 
the 2-D sample of detectable NS-NS binaries for the pa- 
rameters (cost, Di). Since we use an unrealistic mass 
distribution (1.4M -1.4M NS-NS and 1.4M Q -10M Q 
NS-BH binaries), instead of a more astrophysically real- 
istic distribution, the implicit mass prior is uninterest- 
ing. In addition, we assume uniform prior densities for 
the parameter tp between (0,7r). Figure [4] presents the 
average total SNR versus the true Dl of our sample of 
detectable NS-NS and NS-BH binaries for our "full" net- 
work (LIGO, Virgo, AIGO, LCGT). Very few detected 
binaries have SNR above 30 for NS-NS, and above 70 for 
NS-BH. It is interesting to note the different detectable 
ranges in luminosity distance between the two popula- 
tions, with the NS-BH binaries detectable to over twice 
the distance. 

We are also interested in seeing the impact that prior 
knowledge of SHB collimation may have on our abil- 
ity to measure characteristics of these events. We note 
that, to date, there exist only two tentative observations 
which suggest that SHBs may have collimated emission 
(iGrupe et al.ll2006llBvirrows et al.|[200llSoderberg et all 
l2006f k we therefore present results for both moderate 
collimation and for isotropic SHB emission. To obtain 
a sample of beamed SHBs, we assume that the burst 
emission is collimated along the orbital angular momen- 
tum axis, where baryon loading is minimized. Follow- 
ing DHHJ06, we use a distribution for cost = v of 
dP/dv oc exp[-(l - v) 2 /2al], with a v = 0.05. This 



corresponds to a beamed binary population with one 
sigma (68%) of its distribution having an opening jet 
angle within roughly 25°. As Fig. [5] shows, we construct 
a beamed subsample by selecting events from the total 
sample of detected events such that its final distribution 
in inclination angle follows dP/dv. Future joint measure- 
ments of SHBs and GW-driven inspirals should enable us 
to constrain the beaming angles by comparing the mea- 
sured rates of these two populations. 

3.3. Markov Chain Monte Carlo approach 



As mentioned in Sec. 13.11 the principle disadvantage 

of the Bayes estimators #bayes an< ^ ^bayes ^ s ^ nc high 
computational cost of evaluating the multi-dimensional 
integrals which define them, Eqs. (|48|) and (|49|) . To 
get around this difficulty, we use Markov Chain Monte 
Carlo (MCMC) methods to explore the PDFs describ- 
ing the seven parameters {M z , Mi ^l-, cos t, tp, t c , 3> c }. 
MCMC methods are widely used in diverse astrophys- 
ical applications, ranging from high precision cosmology 
(e.g. iDunklev et al] 120091 ISievers et al.ll2009l) to extra- 
solar planet studies (e.g. lFordll2005l I Winn et al1l2007h . 
They have seen an increase in use in GW measure- 
ment and parameter estimation studies in recent years, 
both for the c ase of analysis by the space-based detec- 
tor LISA (e.g-IStroeer et al.ll2006LIWickham et aTl l2006l. 
ICornish fc Porterll2007l . iPorter fc Cornishll200cC and for 
analysis by a net work of ground-bas e d detectors at initial 
sensi tivity (e.g., iRover et all 120071 Ivan der Sluvs et all 
120081) . In this section, we briefly introduce MCMC meth- 
ods, focusi ng in particular on the Metropolis-H astings 
algorithm (jMetropolis et~aTlll953l . lHastingslll970f ). Our 
goal is to describe the workings of our algorithm, rather 
than to review MCMC methods generally; we defe r de- 
tailed proofs to the literature on this s ubject (e.g. iNeall 
[19931 iGilks et~a^lT996l . lMacKavil2003h . Our discussion 
owes much to the helpful bac kground provided by Sec. 
Ill of iChristensen et al.1 (l200l . 

As the "Monte Carlo" part of its name implies, the 
central idea of MCMC is to generate a random sequence 
of parameter states that sample the posterior distribu- 
tion, p{6\s). Let the nth sample in the sequence be 9", 
Then, if one draws a total of N random samples, Eqs. 
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Fig. 3. — The 2-D marginalized prior distribution in luminosity 
distance, Dl, and cosine inclination, cos t. Each point represents a 
detected NS-NS binary for a network comprising all five detectors. 
Notice the bias toward detecting face-on binaries (cost — > ±1); 
these are detected to much larger distances than edge-on (cos t — > 
0). 
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Fig. 4. — Average SNR versus luminosity distance of the total 
detected NS-NS and NS-BH binaries using a network of both LIGO 
detectors, Virgo, AIGO, and LCGT. The left panel shows the total 
detected NS-NS binaries (one point with SNR above 100 is omit- 
ted); the right panel shows the total detected NS-BH binaries (one 
point with SNR above 350 is omitted). Notice the different scales 
of both axes between the two panels. In particular, NS-BH binaries 
are detected to more than twice the distance. 



([48]) and ([49]) can be approximated as simple sample av- 
erages: 
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The key to making this technique work is drawing 
a sequence that represents the posterior PDF. The 
Metropolis-Hastings algorithm, which we now briefly de- 
scribe, ensures that, at least asymptotically, the prob- 
ability of selecting a sample 6" is determined by the 
posterior distribution p(8\s). 

Suppose we are currently at step n in our chain of 
sequences; the zeroth step, 6^°', can be selected arbi- 
trarily. We choose a candidate for the next sequence in 
our chain, 8', from the proposal distribution Q{8'\8^ n '). 
(Note that our choice for this function can have a large 
impact on how effectively the chain we build converges 
to a representation of the posterior PDF; we describe our 
particular choice later in this section.) We then accept 
(or reject) 8' by computing the acceptance probability 



a(8';8 {n) ) 



p(0')<3(0 ( ™ ) ;0') 

P(0(™))Q(0';0(")) 



(52) 



where P(9) is shorthand for the posterior PDF p(0\s). 
We also generate a uniform random deviate u G [0,1]. 
If u < a(6']8 n ), then we accept this step, putting 
0(n+i) = 0/ If u > a(8';8 n ), then this step is re- 
jected, and we put f?(™ +1 ) = 0( n >. Note that sample 
#(«+!) on iy depends upon the previous sample, 8^ n \ 
This is the defining condition of a Markov sequence, and 
explains the "Markov Chain" part of MCMC. The lit- 
erature on this subject demonstrates that the Markov 
chain so constructed explores the full range of the sam- 
ple space spanned by p (8\s). We point the reader to 
IChristensen et all ( 
to lNeall ([19931 ) and 
depth discussion. 

For completeness, we list several standard concerns of 
MCMC methods. First, as already mentioned in passing, 



20041) for a pedag ogical review, and 
Gilks et aLl (1 19961 ) for much more in- 



this algorithm's performance depends on the choice of 
proposal density Q(8'-,8). The efficiency of the method 
increases when the proposal distribution resembles the 
underlying probability distributions, particularly when 
treating strongly degenerate parameters. Second, a cer- 
tain so-called "burn-in" period is needed in order for a 
chain to cquilibriatc to its invariant distribution. This 
period allows the chain to explore all regions of high prob- 
ability in the posterior PDF; it is especially important if 
the chain begins in some parameter region of low proba- 
bility. Finally, a chain will of course necessarily be finite 
in any practical computation. Understanding the con- 
vergence properties of finite chains is central to ensuring 
that the samples we generate are sufficiently large, and 
are effectively independent. Some errors enter because 
of correlations between successive elements of the chain 
and the shot noise. We are especially concerned with es- 
tablishing that the MCMC chain has fully sampled the 
distribution, and is not "stuck" on some parameter island 
of high probability. 

The MCMC algorithm which we use is based 
on a generic versio n of CosmoMC 8 , described in 
iLewis fc Bridld (|2002| ) . We state here only the main fea- 
tures of our code; further d etails of our MCMC m ethod 
are discussed at length in ILewis fc Bridld (|2002t). As 
Eq. ([52]) illustrates, computing a(8'; 8^) requires evalu- 
ating the posterior PDF p(8 | s) for any sample flW . This 
in turn requires that we calculate both the prior distri- 
bution p( >(8), and the likelihood function £tot(s|0) 
[cf. Eq. (|40|) ] for 8 = 8^ n \ We now briefly emphasize key 
aspects of our approach. 

As already mentioned, the total likelihood £tot(s | 8) 
depends upon the individual likelihoods at each detec- 
tor, C a (s a | 8). By Eq. (141]) . we see that each C a (s a \ 8) 
is a function of the detector's output s a (t), the predicted 
waveform h a (6^) at 8^ n \ and the noise spectral density 
S n (f). For a given source (i.e., a given true parameteri- 
zation 0), we need only compute s a {t) once [cf. Eq. 1(53])], 
When computing steps in our chain, we must then com- 




-0.5 0.0 

cos (inc) 

Fig. 5. — Distributions in inclination, cost, for the subsam- 
ple of NS-NS binaries for which we assume SHB collimation 
(dashed line) and for the full sample of total detected binaries 
with isotropic distribution in inclination (solid line). The beamed 
subsample of binaries have a distribution in cos i = v given by 
dP/dv <x exp[-(l - i>) 2 /2o-2], with a v = 0.05. 



8 See |http:7 / cosmologist .info/ cos momc / 1 
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pute the model signal h a (6^) at each sample 6^ n \ 

Priors play a crucial role in our MCMC approach. We 
take the prior distributions in chirp mass M. z , reduced 
mass n z , polarization angle ip, coalescence time t Cl and 
coalescence phase $ c to be flat over the region of sample 
space where the binary is detectable according to our 
selection procedure. More specifically, we choose 

• p^\M z ) = constant over the range [1 M , 2M ] 
for NS-NS; and over the range [2.5 M Q , 4.9 M ] 
for NS-BH. (Note that the true chirp masses in 
the binaries' rest frames are 1.2 M for NS-NS and 
3.0 M Q for NS-BH.) 

• p^°'(fJ>z) = constant over the range [0.3 Mq, 2Mq] 
for NS-NS; and over the range [0.5 Mq, 3.5 M ] 
for NS-BH. (Note that the true reduced masses in 
the binaries' rest frames are 0.7 M for NS-NS and 
1.2 M e for NS-BH.) 

• p(°'(ip) = constant over the range [0, n]. 

• (t c ) — constant over the range 
[—100 sec, 100 sec]. Since we assume that the 
coalescence time is close to the time of the SHB 
event, t c is essentially the time offset between the 
system's final GWs and its SHB photons. We find 
that our choice of prior here is almost irrelevant, 
as long as the prior is fiat and includes the true 
value. No matter how broad we choose the prior 
in t c , our posterior PDF ends up very narrowly 
peaked around t c . 

• p(°)($ c ) = constant over the range [0, 271"]. 

Given our assumed constant comoving density of SHBs, 
the prior distribution for luminosity distance scales as 
comoving volume over the range [0, 2 Gpc] for NS-NS 
binaries, and over the range [0, 5 Gpc] for NS-BH bina- 
ries. For our sample with isotropic inclination distribu- 
tion, we put p(°)(cost) = constant over the range [— 1, 1]. 
For our sample that assumes SHB collimation, our prior 
in cosi = v is the same as the one that we used in our 
selection procedure discussed in the previous subsection: 



dp™ 
dv 



(v) oc e-W** , 



(53) 



with a v = 0.05. 

We then map out full distributions for each of our seven 
parameters, assessing the mean values [Eq. (|4"5|) ] and the 
standard deviations [Eq. (|4"9"|) ]. We generate four chains 
which run in parallel on the CITA "Sunnyvale" Cluster. 
Each chain runs for a maximum of 10 7 steps; we find 
that the mean and median number of steps are ~ 10 5 
and ~ 10 4 , respectively. Each evaluation of the likeli- 
hood function takes ~ 0.3 seconds. We use the first 30% 
of a chain's sample states for "burn in," and thus dis- 
card that data. When generating our Markov chains, we 
use an initial Gaussian proposal distribution with a stan- 
dard deviation similar to that expected for the posterior 
probability density: 



Q{0'-0 



1 



(2ttct t ) w / 2 



exp[-0 2 /2c4], 



(54) 



where o~t is the width of the trial Gaussian distribution 
and N is the number of dimensions of our parameter set 
(in our case, N = 7). However, in the case of parame- 
ters such as Dl and cost that we expect to be strongly 
correlated, we use a much smaller proposal standard de- 
viation. The proposal distribution is then updated using 
the covariance matrix of the last half of the samples. 

Our chains start at random offset parameter values, 
drawn from Gaussians centered on the true parameter 
value. We assess convergence by testing whether the mul- 
tiple chains have produced consistent parameter distribu- 
tions. Following standard practice, we use the Gelman- 
Rubin convergence criterion, defining a sequence as "con- 
verged" if th e statistic R < 1. 1 on t he last half of our 
samples; see iGelman &: Rubinl (|1992l) for more details. 
We use convergence as our stopping criterion. 

Each simulation for every binary runs for an hour to 
forty-eight hours; the mean and median runtime arc eight 
and three hours, respectively. CosmoMC features such as 
temperature annealing and sophisticated sampling meth- 
ods have not been required in our simulations. 

3.4. The "averaged" posterior PDF 

Central to the procedure outlined above is the use of 
the datastream s = h(0) + n which enters the likeli- 
hood function £tot(s|0). The resulting posterior PDF, 
and the parameters one infers, thus depend on the noise 
n which one uses, either via the noise in a particular 
measurement or (as in our case) by simulation (drawing 
randomly from the noise's PDF). In some cases, one will 
want to evaluate statistics that are in a well-defined sense 
"typical" given the average noise properties, rather than 
depending on a particular noise instance. Such averaging 
is appropriate, for example, when forecasting how well an 
instrument should be able to measure the properties of 
a source or process. As we discuss in more detail in the 
following Section, we have also found it is necessary to 
average when trying to compare our MCMC code's out- 
put with previously published work. 

As explained in detail below, the averaged posterior 
PDF takes a remarkably simple form: It is just the 
"usual" posterior PDF, Eq. (|39|) with the noise n set 
to zero. We emphasize that this docs not mean that one 
ignores noise when constructing the averaged PDF. For 
example, one still meaningfully relates the amplitude of 
the signal to the amplitude of typical rms noise in a detec- 
tor by the average SNR, Eq. (|45|) . As such, the averaged 
statistics will show an improvement in measurement ac- 
curacy as the SNR is increased. 

To develop a useful notion of averaged posterior PDF, 
consider the hypothetical (and physically unrealistic) 
case in which we measure a signal using M different noise 
realizations for the same event. The joint posterior PDF 
for these measurements is 

M 

Pjoint(0|si,s 2 , • ■ • sat) = n^ls*) ■ ( 55 ) 

i=l 

Let us define the "average" PDF as the geometric mean 
of the PDFs which describe these measurements: 

p avc (0|s) =p joint (0|s 1 ,s 2 ,...SA/) 1/M . (56) 
Expanding this definition, we find 

M 

p avc (0\s) = Y[lp(8\s t )] 1/M 
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= AAp(°)(0)n^TOT(s,|0)] 1/M , (57) 

i=l 

where the subscript i denotes the ith noise realization in 
our set of M observations. The multi-observation likeli- 
hood can in turn be expanded as 

M M 
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(58) 



x TT cxp 

i=l 

By taking M to be large, the products on the last two 
lines of Eq. (|58|) can be evaluated as follows: 
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h a (6)-h a (0] 



h a (9) - fc o (0) 



= 1 . (59) 

Here, (. . .) denotes an ensemble average over noise real- 
izations (cf. Sec. l3.1[) . and we have used the fact that our 
noise has zero mean. Similarly, we find 
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(60) 

The final equality uses ((n a \n a )) = 2, which can be 
proved using using the noise properties (|3"5)) . (f3"6"|) . and 

(EZD- 

Putting all this together, we finally find 

Pavc (0|s) J] e -(".c)-"-(») I , 

a 

(61) 

where we have absorbed the factor e into the normal- 
ization TV. T/ie posterior PDF, averaged over noise real- 
izations, is simply obtained by evaluating Eq. i39\) with 
the noise n set to zero. 



4. RESULTS I: VALIDATION AND TESTING 

We now validate and test our MCMC code against pre- 
viously published results from CF94. In particular, we 
examine the posterior PDF for a particular NS-NS binary 
which was studied in detail in CF94. We also explore 
the dependence of distance measurement accuracies on 
the detector network and luminosity distance, focusing 
on the strong degeneracy that exists between cos i and 
D L . 

4.1. Comparison with CF94 

Validation of our MCMC results requires comparing 
to work which goes beyond the Gaussian Fisher matrix 
approximation to the likelihood function. In Section IVD 
of CF94, Cutler & Flanagan investigate dominant effects 
that are non- linear in 1/SNR, and consequently, in the 
variable 1/Dt,. Their work shows that including such 
effects has a significant impact on the predicted distance 
measurement accuracies, particularly in the limit of low 
SNR. In particular, they find that Fisher-based estimates 
understate distance measurement errors for a network 
comprising the two LIGO detectors and Virgo. 

Because they go beyond a Fisher matrix analysis, the 
results of CF94 arc a useful comparison to our MCMC 
results. Their paper is also a useful touchstone in that 
they take a binary's source position to be known. It 
should be noted that their motivation for this assump- 
tion is different from ours: while we assume that the 
inspiral is associated with an SHB which determines the 
sky position, they argue that GWs on their own deter- 
mine sky position precisely enough to break correlations 
between position and Dz,. This argument is b ased on 
the "Markovic approximation" ( M arkovid 119931 ) . which 
argues that timing information between different detec- 
tors fixes the sky position to sufficient accuracy. 

Our approach is sufficently different from CF94 that 
we do not expect perfect agreement between our re- 
sults. The most important difference is that we directly 
map out the posterior PDF and compute sample aver- 
ages using Eqs. (|48|) and ([49]) . for all of our parameters 
{M z , Hz, Dl, cos l, ip, t c , $ c }. In contrast, CF94 estimate 
measurement errors for only one parameter, Dl, using 
an approximate method based on an analytical Bayesian 
derivation of the marginalized PDF for Dl ■ Specifically, 
Cutler & Flanagan expand the exponential factor in Eq. 
(|3"9"|) beyond second order in terms of some "best-fit" 
maxiumum likelihood parameters. Their approximation 
treats strong correlations between the parameters Dl 
and cost that are non-linear in 1/SNR. However, other 
correlations between Dl and (tp,(f> c ) are only considered 
up to linear order. They obtain an analytical expression 
for the posterior PDF of the variables Dl and cost in 
terms of their "best-fit" maximum-likelihood values Dl 
and cos! [see Eq. (4.57) of CF94]. The marginalized 1-D 
posterior PDFs for Dl are then computed by numeri- 
cally integrating over the parameter cost. We use the 
term 1-D marginalized PDF in parameter 9i to refer to 
the distribution 
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where p(9\s) is the posterior PDF given by Eq. (f3"9"| and 
N is the number of dimensions of our parameter set (in 
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our case, N = 7). 

In addition to this rather significant difference in our 
techniques, there are several relatively minor differences 
which also affect our comparison: 

• We use the restricted 2PN waveform; CF94 use the 
leading "Newtonian, quadrupole" waveform (cf. 
the waveform we use for pedagogical purposes in 
Sec. II. 2p . Since the distance is encoded in the 
waveform's amplitude, we do not expect that our 
use of a higher-order phase function will have a 
large impact. However, to avoid any easily cir- 
cumvented mismatch, we adopt the Newtonian- 
quadrupole waveform for the purpose of our com- 
parisons. Note that this waveform does not depend 
on the reduced mass /i; hence, for the purpose of 
this comparison only, our parameter space is re- 
duced to 6 dimensions. 

• We use the projected advanced sensitivity noise 
curve shown in Fig. fTJ); CF94 use an analytical 
form [their Eq. (2.1) 9 ] based on the best-guess for 
what advanced sensitivity would achieve at the 
time of their analysis. Compared to the modern 
projected sensitivity, their curve underestimates 
the noise at middle frequencies (~ 40— 150 Hz) and 
overestimates it at high frequencies (> 200 Hz). 
We likewise adopt their noise curve for our com- 
parison. Because of these differences, CF94 rather 
seriously overestimates the SNR for NS-NS inspi- 
ral. Using their noise curve, the average SNR for 
the binary presented in their Fig. 10 is 12. 4 10 ; using 
our up-to-date model for advanced LIGO, it is 5.8. 
As such, the reader should view the numbers in this 
section of our analysis as useful only for validation 
purposes. 

• The two analyses use somewhat different prior dis- 
tributions. As extensively discussed in Sec. 13.31 we 
set uniform priors on the chirp mass M. z , on the 
time t c and phase <f> c at coalescence, and on the 
polarization phase ip. For the purpose of this com- 
parison, we will assume isotropic emission, so we 
also set a flat prior on the inclination cos i. We as- 
sume our sources are uniformly distributed in con- 
stant comoving volume. However, we set a detec- 
tion threshold which depends on the total network 
SNR, which effectively sets a joint prior on source 
inclination and distance. CF94 use a prior distri- 
bution only for the parameter set {Dl, cos l, ijj, <3> c } 
that is flat in polarization phase, coalescence phase, 
and inclination. [This smaller prior set follows from 
their use of the Markovic approximation, which 
also argues that correlations between intrinsic and 
extrinsic parameters can be subsumed into cor- 
relations with the single intrinsic parameter $ c .] 
They assume a prior that is uniform in volume, 
but that cuts off the distribution at a distance 

Anax - 6.5 GpC. 

9 Note that it is missing an overall factor of 1/5 (E. E. Flanagan, 
private communication). 

10 CF94 actually report the SNR to be 12.8. The slight discrep- 
ancy is due to rounding errors in the parameter ro in their Eq. 
(4.28). Adjusting to their preferred value (rather than computing 
ro directly) gives perfect agreement. 



Our goal here is to reproduce the 1-D marginalized 
posterior PDF in Dl for the binary shown in Fig. 10 of 
CF94. For brevity, we refer to this system as the "CF 
binary." Each NS in the CF binary has m z = 1.4 M Q ; it 
has a sky position (6, <f>) = (50°, 276°); and the detector 
network comprises LIGO Hanford, LIGO Livingston and 
Virgo. CF94 report the "best-fit" maximum-likelihood 
values (D L , cosT, *) to be (432 Mpc, 0.31, 101.5°), where 

= tp + Aip(n), and where A^(n) depends on the pre- 
ferred basis of e x and e x set by the detector network 
[see Eqs. (4.23)-(4.25) of CF94 11 ]. In their Bayesian 
framework, the GW signal has been generated and an 
experimenter has no access to the true parameters; they 
thus do not present the parameters' true values. In order 
to compare our distribution with theirs, we assume that 
Q = #ml for the purpose of computing the likelihood 
function C(6\s). This is a reasonable assumption in the 
limit in which priors are taken to be uniform over the rel- 
evant parameter space. As already mentioned, for this 
comparison we use their advanced detector noise curve 
and the Newtonian-quadrupole waveform. Finally, we in- 
terpret the solid curve in Fig. 10 of CF94 as the marginal- 
ized 1-D posterior PDF in Dl for an average of posterior 
PDFs of parameters (given an ensemble of many noisy 
observations for a particular event). We compute the av- 
erage PDF as described in Sec. 13.41 and then marginalize 
over all parameters except Dl, as in Eq. (f62|) . 

The left hand panels of Fig. [6] show the resulting 1- 
D marginalized PDF in Dl and cost . Notice that its 
shape has a broad structure not dissimilar to the solid 
curve shown in Fig. 10 of CF94: The distribution has 
a small bump near Dl ~ 460 Mpc, a main peak at 
Dl ~ 700 Mpc, and extends out to roughly 1 Giga- 
parsec. Because of the broad shape, the Bayesian mean 
(-Dl.bayes = 694 Mpc) is significantly different from 
both the true value (Dl — 432 Mpc in our calculation) 
and from the maximum likelihood (Dl, ml = 495 Mpc). 
Note that, thanks to the marginalization, the peak of this 
curve does not coincide with the maximum likelihood. 
(Interestingly, we find a shape much closer to CF94 Fig. 
10 if we use a prior in which binaries are uniformly dis- 
tributed in distance, rather than uniformly distributed 
in volume.) 

We further determine the 2-D marginalized posterior 
PDFs in Dl and cos l for the CF binary. Fig. [6] illus- 
trates directly the very strong degeneracy between these 
parameters, as expected from the form of Eqs. (l7l) and 
©, as well as from earlier works (e.g., rMarkovi(j ll993l . 
CF94). It's worth noting that, as CF94 comment, the bi- 
nary they chose is measured particularly poorly. This is 
largely due to the fact that one polarization is measured 
by the network far better than the other, so that the 
Dl-cosl degeneracy remains relatively unbroken. This 
degeneracy is responsible for the characteristic tail to 
large Dl we find in the 1-D marginalized posterior PDF 
in Dl, p(Dl\s), which we investigate further in the fol- 
lowing section. 

11 Note that Eq. (4.25) of CF94 should read tan(4A</>) = 
2e +x /(e ++ - e X x)- In addition, § = 56.5° should read 
'I' = 101.5° under the caption of Fig. 10. (We have changed nota- 
tion from ip in CF94 to \E> to avoid multiple accents on the best fit 
value.) We thank Eanna Flanagan for confirming these corrections. 
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Fig. 6. — The 1-D and 2-D marginaliz ed p osterior PDFs for Dl 
and cos t, averaged (as described in Sec. I3.4|l over noise ensembles, 
for the "CF binary." The purpose of this calculation is to repro- 
duce, as closely as possible, the non-Gaussian limit summarized in 
Fig. 10 of CF94. Top left-hand panel shows the 1-D marginalized 
posterior PDF in Dl (the true parameter value Dl = 432 Mpc is 
marked with a solid black line); bottom left-hand panel illustrates 
the 1-D marginalized posterior PDF in cos t (where the true value 
cost = 0.31 is indicated with the solid black line). The right-hand 
panel shows the 2-D marginalized posterior PDF for the parame- 
ters Dl and cos t, where the true values Dl and cos I are marked 
with a cross. The contours around the dark and light blue areas 
indicate the 68 and 95% interval levels, respectively. Notice that 
the true values lie within the parameter interval defined by the 
contours of the 68% region. The Bayesian mean and rms measure- 
ment accuracies are (694.4 Mpc, 0.70) and (162 Mpc, 0.229) for 
{Dl, cost), respectively. 

4.2. Test 1: Varying luminosity distance and number of 

detectors 

We now examine how well we measure as a func- 
tion of both the distance to the CF binary and the 
properties of the GW detector network. Figures [7] 
and [H] show the 1-D and 2-D marginalized posterior 
PDFs in D L and cos t for the CF binary at D L = 
{100, 200, 300, 400, 500, 600} Mpc. For all these cases we 
keep the binary's sky position, inclination, and polariza- 
tion angle fixed as in Sec. 14.11 The average SNRs we 
find for these six cases are (going from D L = 100 Mpc 
to 600 Mpc) 53.6, 26.8, 17.9, 13.4, 10.7, and 8.9 (scaling, 
as expected, as l/D L )- Interestingly, the marginalized 
PDFs for both distance and cost shown in Figs. [7] and [8] 
have fairly Gaussian shapes for D L = 100 and 200 Mpc, 
but have very non-Gaussian shapes for D L > 300 Mpc. 
This can be considered "anecdotal" evidence that the 
Gaussian approximation for the posterior PDF breaks 
down at SNR < 25 or so, at least for this case. For lower 
SNR, the degeneracy between cost and D L becomes so 
severe that the 1-D errors on these parameters become 
quite large. 

Next, we consider how measurement accuracy depends 
on properties of the GW detector network. Figure [9] 
shows the 1-D marginalized posterior PDFs in D L for the 
CF binary with different networks. All "true" parame- 
ters are chosen exactly as described in Sec. 14.11 Adding 
detectors to the network does not substantially increase 
the total SNR; we increase the average total SNR from 



12.4 to 14.6 (adding only AIGO), to 12.4 (adding only 
LCGT; its contribution is so small that the change is 
insignificant to the stated precision), or to 14.7 (adding 
both AIGO and LCGT). This change is not enough to 
counter the D L — cos I degeneracy. While this degeneracy 
remains effective, the distance errors remain large and, in 
this case, biased. We remark also that the CF binary has 
a relatively small SNR for LCGT and Virgo: the average 
SNR in our detectors is 8.23 for LIGO-Hanford, 8.84 for 
LIGO-Livingston, 2.91 for Virgo, 8.71 for AIGO, and 1.1 
for LCGT. This pathology of the CF binary is an exam- 
ple of a fairly general trend that we see. As we show in 
Sec. even considering general binaries, randomly dis- 
tributed on the sky and randomly oriented, we often find 
that SNR is quite low in one or more detectors. 

4.3. Test 2: Varying source inclination 

One of the prime results seen in our analysis of the CF 
binary is (as expected) a strong degeneracy between cos t 
and D L - As Fig. [6] shows, this results in a tail to large 
values of D L in the 1-D marginalized posterior PDF in 
Dl, p(Dl\s), with a Bayes mean Dl for the CF binary 
of Dl = 694 Mpc (compared to the true CF value of 
432 Mpc) . Such a bias is naturally a great concern for 
using these sources as standard sirens, as well as GW 
measurements in general. 

The CF binary has cost ~ 0.31, meaning that it is 
nearly edge-on to the line of sight. Hypothesizing that 
the large tails may be due to its nearly edge-on nature, 
we consider a complementary binary that is nearly face 
on: We fix all of the parameters to those used for the 
CF binary, except for the inclination, which we take to 
be cost = 0.98. We call this test case the "face-on" 
CF binary. By changing the inclination to a more nearly 




Fig. 7. — 1-D and 2-D mar gina lized PDFs for Dl and cos t, 
averaged (as described in Sec. 13.41 1 over noise ensembles for the 
"CF binary" at different values of true luminosity distance Dl'- 
[100 Mpc, 200 Mpc, 300 Mpc] (top to bottom). The true values 
are marked with the solid black line, or a black cross in the 2- 
D case. The Bayesian means and rms measurement errors of the 
luminosity distance are [101.0 Mpc, 212.1 Mpc, 411.2 Mpc] and [3.6 
Mpc, 21.4 Mpc, 110.0 Mpc], respectively. The corresponding Bayes 
mean and rms measurement errors for cos t are [0.317, 0.357, 0.562] 
and [0.033, 0.089, 0.247]. The dark and light contours in the 2-D 
marginalized PDF plots indicate the 68 and 95% interval levels, 
respectively. Notice that the true value always lies within the 68% 
contour region of the 2-D marginalized area at these luminosity 
distances. 
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Fig. 8. — The 1-D and 2-D m argin alized PDFs for Dl and cost, 
averaged (as described in Sec. 13.41 1 over noise ensembles for the 
"CF binary" at different values of the true luminosity distance Dl : 
[400 Mpc, 500 Mpc, 600 Mpc] (top to bottom). The true values are 
marked with the solid black line, or a black cross in the 2-D case. 
The Bayesian means and rms measurement errors of the luminosity 
distance are [627.17 Mpc, 857.3 Mpc, 1068 Mpc] and [148.8 Mpc, 
198.1 Mpc, 262.2 Mpc], respectively. The Bayes means and rms 
measurement errors for cost are [0.686, 0.745, 0.746] and [0.237, 
0.209, 0.218]. The dark and light contours in the 2-D marginalized 
PDF plots indicate the 68 and 95% interval levels, respectively. 
Notice that the true value lies within the 68% contour region of 
the 2-D marginalized area at a luminosity distance of 400 Mpc; 
for larger luminosity distances, the true value lies within the 95% 
contour region. 

face-on situation, we substantially augment the measured 
SNR; the average SNR for the face-on CF binary mea- 
sured by the LIGO/Virgo base network is 24.3 (as com- 
pared to 12.4 for the CF binary). We thus expect some 
improvement in the posterior PDF simply owing to the 
stronger signal. 

Figure [TO] shows the 1-D and 2-D marginalized poste- 
rior PDFs in Dl and cos t. As we might have guessed, 
these distributions arc complementary to those we found 
for the CF binary. In particular, we find that the peak of 
the 1-D marginalized posterior PDF in Dl is shifted to 
lower values in Dl , and the Bayes mean is much closer to 
the true value: Dl = 376.3 Mpc. Notice that the shape 
of the 1-D marginalized posterior PDF in cos t is abruptly 
cut off by the upper bound of the physical prior cos t < 1. 
The Bayes mean for the inclination is cos i = 0.83. 

Just as we varied distance and detector network for the 
CF binary, we now do so for the face-on CF binary. Fig- 
ureslTl1and[T2lshow the 1-D marginalized posterior PDFs 
in Dl and cos t for the facc-on CF binary for Dl = {100, 
200, 300, 400, 500, 600} Mpc. The SNR in these cases 
is 105.0, 52.5, 34.0, 26.2, 20.1, and 17.5, respectively. 
Notice that the 1-D marginalized PDFs for cos t and Dl 
remain highly non-Gaussian, despite the high SNR. Once 
again, this is due to the hard cutoff in the PDF for cos i. 
As long as a significant fraction of the PDF is ex cluded 
by the cutoff (i.e. as long as 1 — cos t < v / ^1XyW T ) ' then 
the TD marginalized TD PDFs for distance and inclina- 
tion will be asymmetric and non-Gaussian. This implies 
that the Gaussian approximation and the Fisher matrix 
method will be unreliable. Similarly, Fig. ll3l shows the 1- 



D marginalized posterior PDFs m Dl for the CF binary 
with different networks. Our results hence confirm this 
predicted degeneracy for face-on binaries irrespective of 
SNR. 

4.4. Discussion of validation and tests 

The main result from our validation tests is that the 
posterior PDFs we find have rather long tails, with strong 
correlations between cost and Dl- Except for cases with 
very high SNR, the 1-D marginalized posterior PDF in 
cost is rather broad. The Bayes mean for cost thus 
typically suggests that a binary is at an intermediate 
inclination — underestimating cos i for nearly face-on bi- 
naries, and overestimating it for nearly edge-on binaries. 
Thanks to the strong cos l-Dl degeneracy, we thus in- 
terpret (at fixed SNR) a nearly face-on binary as be- 
ing closer than its true value, while interpreting a nearly 
edge-on binary as farther than its true value. 

A bias in the measurement oi Dl may be worrying for 
plans to use GW measurements as standard sirens. How- 
ever, thus far we have only considered the measurement 
of single events. The joint constraints from measuring 
an ensemble of events can significantly reduce the bias, 
and improve the overall measurement. To see what hap- 
pens when many events are observed, we consider the 
(physically unrealistic) case in which 80 NS-NS binaries 
are randomly placed on the sky with random orienta- 
tion, but at fixed luminosity distance, Dl = 432 Mpc, 
and measure these using our LIGO/Virgo baseline net- 
work. As we select binaries at a fixed distance on the 
sky, we now assume a uniform distance prior. We assem- 
ble the joint 1-D posterior likelihood for this particular 
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Fig. 9. — T he 1- D marginalized PDF for Dl, averaged (as de- 
scribed in Sec. 13.41 ) over noise ensembles for the "CF binary" for 
different detector networks. The true values of Dl are marked 
with the solid black line. The lower right panel gives the case for 
the CF binary detected by LIGO, Virgo, AIGO, and LCGT. For 
this case, the Bayes mean and rms measurement error in Dl are 
[660.7 Mpc, 153.1 Mpc], respectively. The upper left is LIGO plus 
Virgo with [694.4 Mpc, 162.0 Mpc]; upper right is LIGO, Virgo, 
and AIGO, with [673.7 Mpc, 152.7 Mpc]; and lower left is LIGO, 
Virgo, and LCGT, with [706.1 Mpc, 164.6 Mpc], Notice that the 
addition of detectors has little impact on improving measurements 
in Dl for this binary. 
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Fig. 10.— The 1-D and 2-D marginalized PDFs for D L and 
cos l, averaged (as described in Sec. I3.4|l over noise ensembles for 
the "face-on" CF binary. The Bayesian mean and rms measure- 
ment accuracies are (376.3 Mpc, 0.83) and (51.3 Mpc, 0.12) for 
(Dl, cost), respectively. The top left-hand panel shows the T 
D marginalized posterior PDF in Dl (the true parameter value 
Dl = 432 Mpc is marked with a solid black line); the bottom left- 
hand panel illustrates the 1-D marginalized posterior PDF in cos t 
(where the true value cos I = 0.98 is indicated with the solid black 
line); and the right-hand panel demonstrates the 2-D marginalized 
posterior PDF for parameters and cos (,, where the true values 
Dl and cos t are marked with a cross. The dark and light contours 
indicate the 68 and 95% interval levels, respectively. Notice that 
the true values lie within the parameter interval defined by the 
contours of the 68% region, 
distribution of events: 



80 

Pjoint(-Di) = JJPavc 

i=l 



(D L \ Sl 



(63) 



where pave (-Dl | Si) is the 1-D marginalized posterior PDF 
constructed from p a ve(#|si) for measurement i\ the lat- 
ter quantity having been averaged with respect to noise 
as discussed in Sec. 13.41 Figure [14] shows the joint 1-D 
marginalized posterior PDF for D L . Notice that the bias 
is eliminated; in essence, averaging over many positions 
and orientations washes out the bias. 

Though instructive, this limit is idealized. Real events 
will be distributed in distance, and a sample of eighty 
events is not expected in the near future. In the following 
section, we survey how well distances can be determined 
for realistic ensembles of GW-SHB events. 

5. RESULTS II: SURVEY OF STANDARD SIRENS 

We now examine how well various detector networks 
can measure an ensemble of canonical GW-SHB events. 
We consider NS-NS systems, with each neutron star hav- 
ing mass rriz = (l+z)1.4 M Q ; and NS-BH systems, with 
masses m^ s = (1 + z)1.4M and = (1 + z)10M Q . 
We examine measurement by the four detector networks 
we have discussed (LIGO and Virgo; LIGO, Virgo, and 
AIGO; LIGO, Virgo, and LCGT; and LIGO, Virgo, 
AIGO, and LCGT). Finally, we consider both isotropic 
orientation and a beamed subsample, imagining that the 
gamma-rays from SHBs are beamed along the poles of 
the binary. 

For this study we randomly choose events from our 
sample of detected NS-NS and NS-BH binaries (where 



the selection is detailed in Sec. 13. 2[) . A rough estimate 
of the distance to which we can detect these events can 
be derived as follows. We set a total detector network 
threshold of 7.5, implying a threshold per detector of 
7/v5 = 3.4 for a five detector network. Further averag- 
ing Eq. (|45p over all sky positions and orientations yields 
(DHHJ06) 



a, sky— ave 




(64) 



where the subscript "sky-ave" denotes that is the noise- 
averaged SNR for detector a, averaged over all sky po- 
sitions and orientations. For a single detector threshold 
of 3.4, we find that a five detector network has an aver- 
age range of about 600 Mpc for NS-NS events, and about 
1200 Mpc for NS-BH events. If SHBs are associated with 
face-on binary inspiral, these numbers are increased by 
a factor y^Jl- 1.12 (DHHJ06). 

If a const ant comoving rat e of 10 SHBs Gpc~ 3 yr _1 
is assumed (|Nakar et al.ll2006t ) , we expect approximately 
6 GW-SHB events per year for isotropically orientated 
NS-NS binary progenitors, and 44 SHBs per year for 
isotropically orientated NS-BH binary progenitors. If 
these events are face-on, the factor 1.12 increases the ex- 
pected rate to 8 NS-NS and 57 NS-BH GW-SHB events 
per year. We stress that this is a rough approximation, 
since there are large uncertainties in the SHB event rate 
and redshift distribution. 

In all cases we build our results by constructing the 
posterior distribution for an event given a unique noise 
realisation at each detector. We keep the noise realiza- 
tion, in a given detector and for a specific binary, con- 
stant as we add other detectors. This allows us to make 
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Fig. 11. — The 1-D marg inalized PDFs for Dl and cos t, averaged 
(as described in Sec. 13.41 1 over noise ensembles, for the face-on 
CF binary at different values of the true luminosity distance D^: 
[100 Mpc, 200 Mpc, 300 Mpc] (top to bottom). These values are 
marked with the solid black line. The Bayesian means and rms 
measurement errors of the luminosity distance are [92.4 Mpc, 179.4 
Mpc, 264.5 Mpc] and [6.51 Mpc, 17.3 Mpc, 30.4 Mpc], respectively. 
The counterpart Bayes mean and rms measurement errors for cos t 
arc [0.901, 0.869, 0.851] and [0.068, 0.092, 0.107]. The dark and 
light contours in the 2-D marginalized PDF plots indicate the 68 
and 95% interval levels, respectively. 
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Fig. 12.— The 1-D and 2-D marginalized PDFs for D L and 
cos t, averaged (as described in Sec. I3.4|l over noise ensembles, for 
the face-on CF binary at different values of the true luminosity 
distance D L : [400 Mpc, 500 Mpc, 600 Mpc] (top to bottom). The 
true values are marked with the solid black line. The Bayesian 
means and rms measurement errors for the luminosity distance 
are [352.2 Mpc, 437.8 Mpc, 524.5 Mpc] and [45.0 Mpc, 63.1 Mpc, 
82.1 Mpc], respectively. The counterpart Bayes mean and rms 
measurement errors for cost are [0.842, 0.830, 0.821] and [0.120, 
0.131, 0.142]. The dark and light contours in the 2-D marginalized 
PDF plots indicate the 68 and 95% interval levels, respectively. 

meaningful comparisons between the performance of dif- 
ference detector networks. 

5.1. NS-NS binaries 

We begin with the case in which we have we have de- 
tected two hundred NS-NS binaries, either isotropically 
distributed in inclination angle or from our beamed sub- 
sample, using a network with all five detectors. Figure 
[TBI shows scatter plots of the distance measurement ac- 
curacies for our unbeamed and beamed events, with each 
panel corresponding to a different detector network. The 
distance measurement error is defined as the ratio of the 
rms measurement error with the true value 12 Dj,'- 



AD L V^^ 1 



D L 



(65) 



Y i d lD l j s computed using (|5ip . Although our small sam- 
ple size precludes making definitive quantitative state- 
ments, we emphasize some general trends in by Fig. [15] 
which are particularly relevant to standard siren con- 
straints on cosmological parameters: 

• The unbeamed total sample and the beamed sub- 
sample separate into two rather distinct distribu- 
tions. As anticipated, the beamed subsample im- 
proves measurement errors in D L significantly, by 
greater than a factor of two or more. This is pre- 
dominantly due to the associated beaming prior in- 
cluded in our analysis for these sources. The beam- 
ing prior constrains our inclination angle, cos t, to 

12 Our definition differs from that given in CF94, their Eq. 
(4.62). There, the distance measurement is described as the ra- 
tio of the r ms m easurement error with the Bayes mean. We prefer 
to use Eq. I|65|l as we are interested primarily in the measurement 
error given a binary at its true luminosity distance. 



~ 3%, thereby breaking the strong D l -cosl de- 
generacy. By contrast, when no beaming prior is 
assumed, we find absolute errors of 0.1-0.3 in cost 
for the majority of events; see Fig. [T6j The strong 
-D^-cos i degeneracy then increases our distance er- 
rors. It's worth noting that a significant fraction of 
binaries randomly selected from our sample have 
0.5 < | cos t| < 1. As discussed in Sec. 13.21 this 
is due to the SNR selection criterion — at fixed dis- 
tance, face-on binaries are louder and tend to be 
preferred. 

• When isotropic emission is assumed, we find a 
large scatter in distance measurement errors for all 
events, irrespective of network and true distance; 
we find much less scatter when we assume a beam- 
ing prior. This is illustrated very clearly by the 
upper-right panel of Fig. 1151 In that panel, we 
show the scatter of distance measurement error 
versus true distance for the LIGO, Virgo, AIGO 
detector network, comparing to the Fisher-matrix- 
derived linear scaling trend found in DHHJ06. For 
the unbeamed case (magenta points and line), our 
current results scatter around the linear trend; for 
the beamed case, most events lie fairly close to the 
trend. This demonstrates rather starkly the fail- 
ure of Fisher methods to estimate distance mea- 
surement accuracy, especially when we cannot set 
a beaming prior. 

• Adding detectors to the network considerably in- 
creases the number of detected binaries, but does 
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Fig. 13. — T he 1- D marginalized PDF for Dl, averaged (as 
described in Sec. 13.41 1 over noise ensembles, for the face-on CF 
binary for different detector networks. The true values are 
marked with the solid black line. The lower right-hand panel gives 
the case for the CF binary detected by LIGO, Virgo, AIGO, and 
LCGT. The Bayes mean and rms measurement error for Dl in 
this case is [379.2 Mpc, 48.0 Mpc], respectively. The upper left- 
hand panel is LIGO plus Virgo with [376.3 Mpc, 51.3 Mpc]; upper 
right is LIGO, Virgo, and AIGO, with [378.4 Mpc, 47.9 Mpc]; 
and lower left is LIGO, Virgo, and LCGT, with [380.6 Mpc, 47.6 
Mpc] . As with the nearly edge-on CF binary, we again find that the 
additional detectors don't greatly improve measurements in D^. 
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Fig. 14.— The 1-D marginalized PDFs for D L for eighty NS-NS 
binaries randomly placed on the sky with random orientation, but 
at fixed luminosity distance, Dl = 432 Mpc. We assume measure- 
ment using the LIGO and Virgo detectors. The thin colored lines 
show the normalized 1-D marginalized posterior PDFs for each 
NS-NS bin ary, where each event has been averaged (as described 
in Sec. I3.4|l over all noise ensembles. The thick blue line illustrates 
the normalized joint 1-D PDF for Dl for the eighty events. The 
mean and standard deviation in the mean of the luminosity dis- 
tance for eighty binaries are 468.1 Mpc and 10.8 Mpc, respectively. 
The standard deviation in the luminosity distance for each binary 
is approximately 10.8 Mpc X \/80 fa 96.6 Mpc. 

not significantly improve the accuracy with which 
those binaries are measured. The increase we see in 
the number of detected binaries is particularly sig- 
nificant for GW-SHB standard sirens. For instance, 
an important application is mapping out the pos- 
terior PDF for the Hubble constant Hq. As the 
number of events increases, the resulting joint pos- 
terior PDF in Hq will become increasingly well con- 
strained. Additional detectors also increases the 
distance to which binaries can be detected. This 
can be seen in Fig. [15] for the LIGO and Virgo net- 
work, our detected events extend to Dl ~ 600 Mpc; 
the larger networks all go somewhat beyond this. 
Interestingly, networks which include the AIGO de- 
tector seem to reach somewhat farther out. 

It is perhaps a bit disappointing that increasing the 
number of detectors in our network does not improve 
measurement accuracy. We believe this is due to two 
effects. First, a larger network tends to detect a larger 
number of weaker signals, and thus the additional bina- 
ries are poorly constrained. Second, the principle limita- 
tion to our measurement accuracy is the Dl~cosl degen- 
eracy. A substantial improvement in distance measure- 
ment accuracy would require truly breaking this degen- 
eracy (e.g., by applying the beaming prior). 

5.2. NS-BH binaries 

We now repeat the preceding analysis for NS-BH bi- 
naries. Figure [T7] shows scatter plots of measurement 
accuracies for unbeamed and beamed NS-BH binaries. 
Wc find similar trends as in the NS-NS case: 



• The unbeamed and beamed samples separate into 
two distinct distributions. Notice, however, that 
outliers exist in measurement errors at high D l for 
several beamed events for all networks. This is not 
too surprising, given we expect beamed sources at 
higher luminosity distances and lower SNR. Such 
events are more likely to deviate from the linear 
relationship predicted by the Fisher matrix. 

• We see substantial scatter in distance measure- 
ment, particularly when isotropic emission is as- 
sumed. As with NS-NS, the scatter is not as se- 
vere when we assume beaming, and in that case 
lies fairly close to a linear trend, as weould be pre- 
dicted by a Fisher matrix. (Note that this trend is 
shallower in slope than for NS-NS binaries, thanks 
to the larger mass of the system.) 

• We do not see substantial improvement in distance 
measurement as we increase the detector network. 
As with NS-NS, adding detectors does increase the 
range of the network; AIGO appears to particu- 
larly add events at large Dl (for both the isotropic 
and beamed samples). However, adding detectors 
does not break the fundamental Dl~cos l degener- 
acy. From our full posterior PDFs, we find absolute 
errors of 0.1-0.3 in cost, very much like in the NS- 
NS case. 



6. SUMMARY DISCUSSION 

In this analysis, we systematically study how well grav- 
itational waves can be used to measure luminosity dis- 
tances, under the assumption that short-hard gamma 
ray bursts are accompanied by binary inspiral. We ex- 
amine two plausible compact binary SHB progenitors, 
and a variety of plausible detector network configura- 
tions. We emphasize that we assume sky position is 
known. We build on the previous study of DHHJ06, 
which used the so-called Gaussian approximation of the 
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Fig. 15. — Distance measurement errors versus the true lu- 
minosity distance for our sample of NS-NS binaries. Colored 
points assume isotropic emission in our priors; black crosses as- 
sume our beaming prior. The lower right-hand panel shows the 
result for our "full" network (LIGO, Virgo, AIGO, and LCGT); 
in this case, 200 unbeamed and 200 beamed binaries are in- 
cluded. Upper left shows LIGO plus Virgo (84 unbeamed and 
90 beamed); upper right is LIGO, Virgo, and AIGO (140 un- 
beamed and 138 beamed); lower left is LIGO, Virgo, and LCGT 
(141 unbeamed and 141 beamed). In the LIGO, Virgo, AIGO 
panel we also show the Fisher-matrix-dcrived linear scaling given 
in DHHJ06: ADl/Dl — 6^/(4.4 Gpc) assuming beaming, and 
ADl/Dl — Sl/(1.7 Gpc) for isotropic emission. 
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Fig. 16. — Inclination angle measurement errors versus true incli- 
nation angle for NS-NS binaries, assuming an isotropic orientation 
distribution. The lower right-hand panel is for the "full" LIGO, 
Virgo, AIGO, LCGT network and includes 200 binaries. Upper 
left is LIGO plus Virgo (81 binaries); upper right is LIGO, Virgo, 
and AIGO (140 binaries); and lower left is LIGO, Virgo, and LCGT 
(141 binaries). 

posterior PDF. This approximation is known to work 
well for large SNR; the limits of its validity, however, are 
not well understood. Since the SNR of events measured 
by ground-based detectors is likely to be of order 10, we 
are concerned that the Gaussian limit may not be valid. 
We examine the posterior PDF for the parameters of 
measured events using Markov Chain Monte Carlo tech- 
niques, which do not rely on this approximation. We also 
introduce a well-defined posterior PDF that is averaged 
over all noise realizations and does not depend solely on 
a particular noise instance. Such a quantity is especially 
useful when we wish to predict how well a detector should 
be able to measure the properties of a source or process. 

We find that the Gaussian description of the likelihood 
substantially underestimates distance measurement er- 
rors. We also find that the main limitation for individual 
standard siren measurements is the strong degeneracy 
between distance to the binary and its inclination to the 
line of sight. Adding detectors to a network does not sub- 
stantially improve the distance measurement for a given 
single event. When we assume that the SHB is isotropic 
(so that we cannot infer anything about the source's incli- 
nation given the GRB measurement), we find that Fisher 
matrix estimates of distance errors are very inaccurate. 
Our distributions show very large scatter about Fisher- 
based predictions. 

The situation is improved dramatically if we can as- 
sume that SHBs are collimated, thereby giving us a prior 
on the orientation of the progenitor binary. By assum- 
ing that the GRBs are preferentially emitted into an 
opening angle of roughly 25°, we find that the distance- 
inclination correlation is strongly broken. In this case, 
the Fisher matrix estimates are much more reasonable, 
giving a fairly good sense of the trend with which dis- 
tances are determined (albeit with a moderate amount 
of scatter about that trend). This illustrates the impor- 
tance of incorporating prior knowledge of our parameters 
into our measurement. 

Our distance measurement results are summarized by 
Fig. H5] (for NS-NS SHB progenitors) and Fig. [17] (for NS- 
BH). Assuming isotropy, we find the distance to NS-NS 
binaries is measured with a fractional error of roughly 
20-60%, with most events in our distribution clustered 
near 20-30%. Beaming improves this by roughly a factor 
of two, and elminates most of the high error events from 
our sample. Similar results are found for NS-BH events, 



with perhaps an improvement of 10% or so. 

At first blush these findings appear disheartening with 
regards to standard siren measurements. We emphasize, 
however, that these results describe the outcome of in- 
dividual siren measurements. When these measurements 
are used as cosmological probes, we will be interested in 
constructing the joint distribution, following observation 
of TV GW-SHB events. As Fig. [H] illustrates, the joint 
distribution found by combining many individual mea- 
surements becomes increasingly well constrained, wash- 
ing out the scatter we see on an event- by- event basis. 
Indeed, preliminary studies show that our ability to con- 
strain Hq improves quite sharply as the number of mea- 
sured binaries is increased. In our most pessimistic sce- 
nario (the SHB is assumed to be a NS-NS binary, with 
no prior on inclination, and measured by the baseline 
LIGO- Virgo network), we find that Hq can be measured 
with ~ 13% fractional error with N = 4, improving to 
~ 5% for N = 15. We stress that this measurement is 
based on absolute determinations of distance using stan- 
dard sirens. (Details of this analysis will be presented in 
a companion paper, currently in preparation.) 

Increasing the number of measured events will thus 
be crucial for getting cosmologically interesting measure- 
ments. To this end, it is important to note that increas- 
ing the number of detectors in our network enables a 
considerable increase in the number of detected binaries. 
This is due to both improvement in sky coverage, and 
in the total detection volume. Going from a network 
which includes all four detectors (LIGO, Virgo, AIGO, 
and LCGT) to our baseline network of just LIGO and 
Virgo entails a ~ 50% reduction in the number of de- 
tected binaries. Eliminating just one of the proposed 
detectors (AIGO or LCGT) leaves us with - 75% of 
the original detected sample. Interestingly, we find that 
AIGO shows a marked improvement in certain areas of 
the sky, resulting in a smaller scatter of measurement 
accuracies for an ensemble of events at low SNR. 

Assumin g an event rate de nsity of 10 SHBs per year 
per Gpc 3 dNakar et al.l 120061 ) , we expect to measure 6 
NS-NS events per year or 44 NS-BH events per year us- 
ing a network with all five detectors. For the LIGO- Virgo 
network, we expect half this rate; for the LIGO-Virgo- 
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Fig. 17. — Distance measurement errors versus the true lumi- 
nosity distance for our sample of NS-BH binaries. Colored points 
assume isotropic emission; black crosses use our beaming prior. 
The lower right-hand panel shows the sample detected by our "full" 
network (LIGO, Virgo, AIGO, LCGT) and includes 250 unbeamed 
and 200 beamed binaries. Upper left is LIGO plus Virgo (117 un- 
beamed and 98 beamed); upper right is LIGO, Virgo, and AIGO 
(180 unbeamed and 147 beamed); and lower left is LIGO, Virgo, 
and LCGT (179 unbeamed and 144 beamed). 
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AIGO or LIGO-Virgo-LCGT network, we expect 3/4 of 
this rate. If SHB collimation can be assumed, the rate 
is further augmented by a factor of 1.12. At this rate, 
we find that one year of observation should be enough 
to measure Hq to an accuracy of ~ 1% if SHBs are dom- 
inated by beamed NS-BH binaries using the "full" net- 
work of LIGO, Virgo, AIGO, and LCGT— admittedly, 
our most optimistic scenario. A general trend we see is 
a network of five detectors (as opposed to our baseline 
LIGO- Virgo network of three detectors) increases mea- 
surement accuracy in Hq by a factor of one and a half; 
assuming that the SHB progenitor is a NS-BH binary 
improves measurement accuracies by a factor of four or 
greater. Errors in H$ are seen to improve by a factor of 
at least two when we assume SHB collimation. 

Aside from exploring the cosmological consequences of 
these results, several other issues merit careful future 
analysis. One general result we found is the importance 
that prior distributions have on our final posterior PDF. 
We plan to examine this in some detail, checking which 
parameters particularly influence our final result, and as- 
certaining what uncertainties can be ascribed to our in- 
ability to set priors on these parameters. It may be pos- 
sible to mitigate the influence of the £>l~cos l degeneracy 
by setting a distance prior that requires our inferred dis- 
tance to be consistent with the SHB's observed redshift. 

Another important issue is that of systematic errors 
in binary modeling. We have used the second-post- 
Newtonian description of a binary's GWs in our analy- 
sis; and, we have ignored all but the leading quadrupole 
harmonic of the waves (the so-called "restricted" post- 
Newtonian waveform). Our suspicion is that a more 
complete post-Newtonian description of the phase would 
have little impact on our results, since such effects are 
not likely to have an impact on the all- important Dl— 
cos i degeneracy. In principle, including additional (non- 
quadrupolc) harmonics could have an impact on this de- 
generacy, since these other harmonics encode different 
information about the inclination angle u. In practice, 
we expect that they won't have much effect on GW-SHB 
measurements, since these harmonics are measured with 
very low SNR (the strongest harmonic is roughly a fac- 
tor of 10 smaller in amplitude than the quadrupole). It 
shouldn't be too difficult to test this, however; given how 
important this degeneracy has proven to be, it could be 
a worthwhile exercise. 
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As discussed previously, we confine our analysis to the 
inspiral part of the waveform. Inspiral waves are ter- 
minated at the presumed innermost stable circular or- 
bit frequency, /isco = (6 3/2 ttM z ). For NS-NS binaries, 
/isco — 1600 Hz. At this frequency, detectors have fairly 
poor sensitivity, and we are thus confident that termi- 
nating the waves has little impact on our results for NS- 
NS systems. However, for our assumed NS-BH binaries, 
/isco — 400 Hz. Detectors have rather good sensitivity 
in this band, so it may be quite important to improve 
our model for the waves' termination in this case. 

Perhaps the most important follow-up would be to in- 
clude the impact of spin. Although the impact of neutron 
star spin is likely to be small, it may not be negligible; 
and, for NS-BH systems, the impact of the black hole's 
spin could be significant. Spin induces precessions in 
the binary which can make the orientation of the orbit, 
L, dynamical. That in turn makes the observed incli- 
nation dynamical, which can break the D^-cosl degen- 
eracy. Van der Sluys et al. (2008) have already shown 
that spin precession physics vastly improves the ability 
of ground-based detectors to determine a source's posi- 
tion on the sky; we arc confident that a similar analysis 
which assumes sky position will find that measurements 
of source distance and inclination can likewise be im- 
proved. 
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